text:
export load, save

"""
        `load(fitsfile::String, IntensityMap)`

This loads in a fits file that is more robust to the various imaging algorithms
in the EHT, i.e. is works with clean, smili, eht-imaging.
The function returns an tuple with an intensitymap and a second named tuple with ancillary
information about the image, like the source name, location, mjd, and radio frequency.
"""
function fileio_load(file::File{format"FITS"}, T::Type{IntensityMap})
    return _load_fits(file.filename, T)
end

function _load_fits(fname, ::Type{IntensityMap})
    img, head = FITS(fname, "r") do f
        if ndims(f[1])==2
            image = copy(read(f[1])')[:,end:-1:1]
        else
            @warn "Currently only stokes I is loaded"
            image = Matrix{Float64}(read(f[1])[:,:,1,1]')
        end

        header = read_header(f[1])
        nx = Int(header["NAXIS1"])
        ny = Int(header["NAXIS2"])
        psizex = abs(float(header["CDELT1"]))*π/180
        psizey = abs(float(header["CDELT2"]))*π/180

        ra = float(header["OBSRA"])
        dec = float(header["OBSDEC"])
        #Get frequency
        freq = 0.0
        if haskey(header, "FREQ")
            freq = parse(Float64, string(header["FREQ"]))
        elseif "CRVAL3" in keys(header)
            freq = float(header["CRVAL3"])
        end

        mjd = 0.0
        if haskey(header, "MJD")
            mjd = parse(Float64, string(header["MJD"]))
        end

        source = "NA"
        if haskey(header,"OBJECT")
            source = string(header["OBJECT"])
        end

        bmaj = 1.0 #Nominal values
        bmin = 1.0
        if haskey(header, "BUNIT")
            if header["BUNIT"] == "JY/BEAM"
                @info "Converting Jy/Beam => Jy/pixel"
                bmaj = header["BMAJ"]*π/180
                bmin = header["BMIN"]*π/180
                beamarea = (2.0*π*bmaj*bmin)/(8*log(2))
                image .= image.*(psizex*psizey/beamarea)
            end
        end

        imap = IntensityMap(image, psizex*nx, psizey*ny)
        info = (source=source, RA=ra, DEC=dec, mjd=mjd, ν=freq)
        return imap, info
    end
end

"""
    `$(SIGNATURES)`
Saves an image to a fits file. You can optionally pass an EHTObservation so that ancillary information
will be added.
"""
function fileio_save(fname::File{format"FITS"}, img::IntensityMap, obs = nothing)
    head = make_header(obs)
    _save_fits(fname.filename, img, head)
end

function make_header(obs)
    if isnothing(obs)
        return (source="NA", RA=0.0, DEC=0.0, mjd=0, freq=0.0)
    else
        return (source=String(obs.source), RA=obs.ra, DEC=obs.dec, mjd=obs.mjd, freq=obs.frequency)
    end
end

function _save_fits(fname::String, image::IntensityMap, head)
    headerkeys = ["SIMPLE",
                  "BITPIX",
                  "NAXIS",
                  "NAXIS1",
                  "NAXIS2",
                  "OBJECT",
                  "CTYPE1",
                  "CTYPE2",
                  "CDELT1",
                  "CDELT2",
                  "OBSRA",
                  "OBSDEC",
                  "FREQ",
                  "CRPIX1",
                  "CRPIX2",
                  "MJD",
                  "TELESCOP",
                  "BUNIT",
                  "STOKES"]
    values = [true,
              -64,
              2,
              size(image, 1),
              size(image, 2),
              head.source,
              "RA---SIN",
              "DEC---SIN",
              rad2deg(image.psizex),
              rad2deg(image.psizey),
              head.RA,
              head.DEC,
              head.freq,
              size(image,1)/2+0.5,
              size(image,2)/2+0.5,
              head.mjd,
              "VLBI",
              "JY/PIXEL",
              "STOKES"]
    comments = ["conforms to FITS standard",
                "array data type",
                "number of array dimensions",
                "",
                "",
                "",
                "",
                "",
                "",
                "",
                "",
                "",
                "",
                "",
                "",
                "",
                "",
                "",
                ""]
    FITS(fname, "w") do hdu
        hdeheader = FITSHeader(headerkeys, values, comments)
        img = copy(image[:,end:-1:1]')
        write(hdu, img, header=hdeheader)
    end
end

"""
$(SIGNATURES)
Load a ThemisPy style ascii EHT observation file.
"""
function load_tpy(file)
    data = readdlm(file, skipstart=1)
    bs1 = Symbol.(getindex.(data[:,5], Ref(1:2)))
    bs2 = Symbol.(getindex.(data[:,5], Ref(3:4)))
    baselines = tuple.(bs1, bs2)
    edata = StructArray{EHTVisibilityDatum{Float64}}(
                visr=float.(data[:,8]),
                visi=float.(data[:,10]),
                error=float.(data[:,9]),
                u=float.(data[:,6])*1e6,
                v=float.(data[:,7])*1e6,
                time=float.(data[:,4]),
                frequency=fill(227e9, size(data,1)),
                bandwidth=fill(4e6, size(data,1)),
                baseline=baselines
            )
    mjd =  Int(modified_julian(UTCEpoch(Int(data[1,2]), 1,1,0)).Δt+float(data[1,3]))
    return EHTObservation(data=edata, mjd=mjd, ra=180.0, dec=0.0, source=Symbol(data[1,1]))
end


text:
export AbstractOEDParam, getProblemSize, getM, getdpMx, getdpMTx, getD, getdpD, getNumberOfDesignParameters

"""
	abstract OED.AbstractOEDParam

	generic type for OED problems. For an example see OEDTypeAParam
"""
abstract type AbstractOEDParam end

"""
m,n,nrhs = getProblemSize(pFor::AbstractOEDParam)

"""
getProblemSize(pFor::AbstractOEDParam) = error("nyi for pFor of type $(typeof(pFor))")

"""
function OED.getM

constructs matrix M(p) for OED problems

Input:

p     - design parameters
pFor  <: AbstractOEDParam (describing the construction)

Output:

M     - forward operator
"""
getM(p::Vector{Float64},pFor::AbstractOEDParam)       = error("nyi for pFor of type $(typeof(pFor))")

"""
function OED.getdpMx

constructs gradient matrix  d_p (M(p)*x)

Input:

p     - design parameters
x     - projection
pFor  <: AbstractOEDParam (describing the construction)

Output:

M     - forward operator
"""
getdpMx(p::Vector{Float64},x::Vector{Float64},pFor::AbstractOEDParam)  = error("nyi for pFor of type $(typeof(pFor))")

"""
function OED.getdpMTx

constructs gradient matrix  d_p (M(p)'*x)

Input:

p     - design parameters
x     - projection
pFor  <: AbstractOEDParam (describing the construction)

Output:

M     - forward operator
"""
getdpMTx(p::Vector{Float64},x::Vector{Float64},pFor::AbstractOEDParam) = error("nyi for pFor of type $(typeof(pFor))")

"""
function OED.getD

get data for current design

Input:

p     - design parameters
pFor  <: AbstractOEDParam (describing the construction)

Output:

M     - forward operator
"""
getD(p::Vector{Float64},pFor::AbstractOEDParam)       = error("nyi for pFor of type $(typeof(pFor))")

"""
function OED.getdpD

gets gradient matrix d_p (D(p)[:,k])

Input:

p     - design parameters
k     - index of data
pFor  <: AbstractOEDParam (describing the construction)

Output:

M     - forward operator
"""
getdpD(p::Vector{Float64},k::Int,pFor::AbstractOEDParam)     = error("nyi for pFor of type $(typeof(pFor))")

"""
function OED.getNumberOfDesignParameters(pFor)
"""
getNumberOfDesignParameters(pFor::AbstractOEDParam)     = error("nyi for pFor of type $(typeof(pFor))")

ext: jl


text:
joinn(v) = join(v, "\n\n")
parenwrap(str1::String, str2::String) = "($str1\n$str2)"

hasname(item, name) = isa(item, GlobalRef) && item.name == name || item == name
hasname(item, names::Array) = any([hasname(item, name) for name in names])
hasname(item, names::Base.KeySet) = any([hasname(item, name) for name in names])

isnothing(item) = isa(item,Nothing) || isa(item,Array) && length(item)==1 && isa(item[1], Nothing)


text:
using LAR

datafile = readcsv("test/svg/test2.lines");
#datafile = readcsv("test/svg/provaWall.lines");
V = reshape(datafile',(size(datafile',1)÷2,size(datafile',2)*2))
len = length(datafile)
EV = collect(reshape(1:(len÷2), 2,(len÷4)))
larview(V,EV)
lineArray = (V,EV)
W,EW = lines2lar(lineArray)
viewexploded(W,EW)
#VV = vertices2vertices(W,EW)
viewLarIndices(W,EW,0.75)
#viewLarIndices(W,EW,0.025)
chains = boundary(W,EW)
operator = boundaryOp(EW,chains)
#println(full(operator))
FW = [sort(collect(Set(vcat([EW[:,abs(e)] for e in face]...)))) for face in chains]
viewLarIndices(W,EW,FW,0.75)

boxes = lar2boxes(W,FW)
parts = boxBucketing(boxes)

V,FV,EV = larFromLines(datafile)

chainAreas(V,EV,chains)


text:
using Pkg, RomeoApp
Pkg.test("RomeoApp")
@show "test finished"
phasefile = abspath(joinpath(@__DIR__, "data", "small", "Phase.nii"))
magfile = abspath(joinpath(@__DIR__, "data", "small", "Mag.nii"))

function test_romeo(args)
    file = tempname()
    args = [args..., "-o", file]
    unwrapping_main(args)
end

args = [phasefile, "-B", "-t", "1:3"]
test_romeo(args)

args = [phasefile, "-m", magfile, "-t", "1:3"]
test_romeo(args)


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: iotsecuretunneling
using AWS.Compat
using AWS.UUIDs

"""
    CloseTunnel()

Closes a tunnel identified by the unique tunnel id. When a CloseTunnel request is received, we close the WebSocket connections between the client and proxy server so no data can be transmitted.

# Required Parameters
- `tunnelId`: The ID of the tunnel to close.

# Optional Parameters
- `delete`: When set to true, AWS IoT Secure Tunneling deletes the tunnel data immediately.
"""
close_tunnel(tunnelId; aws_config::AWSConfig=global_aws_config()) = iotsecuretunneling("CloseTunnel", Dict{String, Any}("tunnelId"=>tunnelId); aws_config=aws_config)
close_tunnel(tunnelId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iotsecuretunneling("CloseTunnel", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tunnelId"=>tunnelId), args)); aws_config=aws_config)

"""
    DescribeTunnel()

Gets information about a tunnel identified by the unique tunnel id.

# Required Parameters
- `tunnelId`: The tunnel to describe.

"""
describe_tunnel(tunnelId; aws_config::AWSConfig=global_aws_config()) = iotsecuretunneling("DescribeTunnel", Dict{String, Any}("tunnelId"=>tunnelId); aws_config=aws_config)
describe_tunnel(tunnelId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iotsecuretunneling("DescribeTunnel", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tunnelId"=>tunnelId), args)); aws_config=aws_config)

"""
    ListTagsForResource()

Lists the tags for the specified resource.

# Required Parameters
- `resourceArn`: The resource ARN.

"""
list_tags_for_resource(resourceArn; aws_config::AWSConfig=global_aws_config()) = iotsecuretunneling("ListTagsForResource", Dict{String, Any}("resourceArn"=>resourceArn); aws_config=aws_config)
list_tags_for_resource(resourceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iotsecuretunneling("ListTagsForResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceArn"=>resourceArn), args)); aws_config=aws_config)

"""
    ListTunnels()

List all tunnels for an AWS account. Tunnels are listed by creation time in descending order, newer tunnels will be listed before older tunnels.

# Optional Parameters
- `maxResults`: The maximum number of results to return at once.
- `nextToken`: A token to retrieve the next set of results.
- `thingName`: The name of the IoT thing associated with the destination device.
"""
list_tunnels(; aws_config::AWSConfig=global_aws_config()) = iotsecuretunneling("ListTunnels"; aws_config=aws_config)
list_tunnels(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iotsecuretunneling("ListTunnels", args; aws_config=aws_config)

"""
    OpenTunnel()

Creates a new tunnel, and returns two client access tokens for clients to use to connect to the AWS IoT Secure Tunneling proxy server.

# Optional Parameters
- `description`: A short text description of the tunnel.
- `destinationConfig`: The destination configuration for the OpenTunnel request.
- `tags`: A collection of tag metadata.
- `timeoutConfig`: Timeout configuration for a tunnel.
"""
open_tunnel(; aws_config::AWSConfig=global_aws_config()) = iotsecuretunneling("OpenTunnel"; aws_config=aws_config)
open_tunnel(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iotsecuretunneling("OpenTunnel", args; aws_config=aws_config)

"""
    TagResource()

A resource tag.

# Required Parameters
- `resourceArn`: The ARN of the resource.
- `tags`: The tags for the resource.

"""
tag_resource(resourceArn, tags; aws_config::AWSConfig=global_aws_config()) = iotsecuretunneling("TagResource", Dict{String, Any}("resourceArn"=>resourceArn, "tags"=>tags); aws_config=aws_config)
tag_resource(resourceArn, tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iotsecuretunneling("TagResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceArn"=>resourceArn, "tags"=>tags), args)); aws_config=aws_config)

"""
    UntagResource()

Removes a tag from a resource.

# Required Parameters
- `resourceArn`: The resource ARN.
- `tagKeys`: The keys of the tags to remove.

"""
untag_resource(resourceArn, tagKeys; aws_config::AWSConfig=global_aws_config()) = iotsecuretunneling("UntagResource", Dict{String, Any}("resourceArn"=>resourceArn, "tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iotsecuretunneling("UntagResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceArn"=>resourceArn, "tagKeys"=>tagKeys), args)); aws_config=aws_config)


text:
using AstroUtils, StaticArrays

# Test 1: Eliptical non-equatorial

# Position and velocity
rv = SVector(6524.834, 6862.875,  6448.296, # [km]
             4.901327, 5.533756, -1.976341) # [km/s]

# Gravitational parameter
μ = 398600.4418 # [km³/s²]

# Convert
(kep,flag) = AstroUtils.cartToKep(rv, μ)

# Tests
@test abs(kep[1] - 36127.3376) < 0.01
@test abs(kep[2] - 0.8328534)  < 0.0001
@test abs(kep[3] - 1.5336056)  < 0.001
@test abs(kep[4] - 3.9775750)  < 0.001
@test abs(kep[5] - 0.9317428)  < 0.0001
@test abs(kep[6] - 1.6115525)  < 0.001
@test flag == 0


text:
module Plane

function get_fn(argd::Dict)
    center = [float(x) for x in argd["center"]]
    raw_direction = [float(x) for x in argd["direction"]]
    half_height = 0.5 * float(argd["height"])
    k = float(argd["k"])
    power = convert(Int, argd["power"])

    direction = raw_direction/norm(raw_direction)

    function calc(d::AbstractFloat, doCalc::Bool)
        if doCalc
            return (1.0/power) * k * (d^power)
        else
            return 0.0
        end
    end

    function fn_compute(coordinates)
        d = abs(dot(coordinates - center, direction))
        return calc(d, d > half_height)
    end

    return fn_compute
end

export get_fn
end

max_stars_repo_name: yuhangwang/SimShape

max_issues_repo_name: yuhangwang/SimShape

max_forks_repo_name: yuhangwang/SimShape


text:
module Issue722

using Test
using Gridap
using Gridap.Fields
using Gridap.CellData

model = simplexify(CartesianDiscreteModel((0,1,0,1),(2,2)))
Ω=Triangulation(model)
Ωo=Gridap.Geometry.TriangulationView(Ω,[1,2,3])
glue = get_glue(Ωo,Val(2))
@test glue.tface_to_mface == [1, 2, 3]
@test glue.mface_to_tface == [1, 2, 3, -1, -2, -3, -4, -5]
no=GenericCellField(fill(GenericField(x->x[2]),num_cells(Ωo)),Ωo,ReferenceDomain())
x_Ω = get_cell_points(Ω)
nox = no(x_Ω)
collect(nox)

end # module


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)
# The above causes the warning: eval from module __anon__ to Proj4

module Proj4
using Compat
import Compat.String

const _projdeps = joinpath(dirname(@__FILE__),"..","deps","deps.jl")
if isfile(_projdeps)
    include(_projdeps)
else
    error("Proj4 is not properly installed. Please run Pkg.build(\"Proj4\")")
end

export Projection, # proj_types.jl
       transform, transform!,  # proj_functions.jl
       is_latlong, is_geocent, compare_datums, spheroid_params,
       xy2lonlat, xy2lonlat!, lonlat2xy, lonlat2xy!

include("projection_codes.jl") # ESRI and EPSG projection strings
include("proj_capi.jl") # low-level C-facing functions (corresponding to src/proj_api.h)

function _version()
    m = match(r"(\d+).(\d+).(\d+),.+", _get_release())
    VersionNumber(parse(Int, m[1]), parse(Int, m[2]), parse(Int, m[3]))
end

"Parsed version number for the underlying version of libproj"
const version = _version()

# Detect underlying libproj support for geodesic calculations
const has_geodesic_support = version >= v"4.9.0"

if has_geodesic_support
    export geod_direct, geod_inverse, geod_destination, geod_distance
    include("proj_geodesic.jl") # low-level C-facing functions (corresponding to src/geodesic.h)
end

include("proj_types.jl") # type definitions for proj objects
include("proj_functions.jl") # user-facing proj functions

"Get a global error string in human readable form"
error_message() = _strerrno()

end # module


text:
###################################################
#
# Copyright 2016 SOASTA, Inc.
# Distributed under the terms of the MIT license
#
# File: RepositoryAPI.jl
#
# Functions to communicate with the mPulse Repository REST API
# This file MUST be `include()`d from `mPulseAPI.jl`
#
###################################################

export
    deleteRepositoryObject,
    getRepositoryObject,
    postRepositoryObject

const TokenTimeoutHours = 5

# Internal convenience function.  Fetches an object from the repository and caches it for an hour in the appropriate cache object
# - Returns a single object if filter keys are passed in an filterRequired is set to true (default)
# - Returns an array of objects if filter keys are not passed in and filterRequired is set to false
# - Throws an exception if filter keys are not passed in and filterRequired is set to true
function getRepositoryObject(token::AbstractString, objectType::AbstractString, searchKey::Dict{Symbol, Any}; filterRequired::Bool=true)
    global verbose

    if token == ""
        throw(ArgumentError("`token' cannot be empty"))
    end

    local isKeySet = false

    if filterRequired
        isKeySet = any(kv -> isa(kv[2], Number) ? kv[2] > 0 : !isempty(kv[2]), searchKey)

        if !isKeySet
            throw(ArgumentError("At least one of `$(join(collect(keys(searchKey)), "', `", "' or `"))' must be set"))
        end
    end

    local object = getObjectFromCache(objectType, searchKey)

    if object != nothing
        return object
    else
        object_list = getHttpRequest(token, objectType, searchKey, isKeySet)
    end

    if isKeySet
        return object_list[1]
    else
        return object_list
    end
end

# Internal convenience function.  Updates an object from the repository.
function postRepositoryObject(token::AbstractString,
                              objectType::AbstractString,
                              searchKey::Dict{Symbol, Any};
                              attributes::Dict=Dict(),
                              objectFields::Dict=Dict(),
                              body::Union{AbstractString, LightXML.XMLElement}="",
                              filterRequired::Bool=true
)

    global verbose

    if token == ""
        throw(ArgumentError("`token' cannot be empty"))
    end

    objectID = get(searchKey, :id, 0)
    name = get(searchKey, :name, "")

    object = getObjectInfo(token, objectType, objectID, name)

    # If objectID was not supplied, it will now be available
    objectID = get(object, "id", 0)

    # Retrieve existing (old) attributes
    oldAttributes = Dict{AbstractString, Any}(get(object, "attributes", Dict()))

    local isKeySet = false

    if filterRequired
        isKeySet = any(kv -> isa(kv[2], Number) ? kv[2] > 0 : !isempty(kv[2]), searchKey)

        if !isKeySet
            throw(ArgumentError("At least one of `$(join(collect(keys(searchKey)), "', `", "' or `"))' must be set"))
        end
    end

    local url = ObjectEndpoint * "/" * objectType * "/$(objectID)"
    local debugID = "(all)"

    if verbose
        println("POST $url")
        println("X-Auth-Token: $token")
        !isempty(attributes) && println(attributes)
        !isempty(objectFields) && println(objectFields)
    end

    json = buildPostJSON(objectType, objectID, objectFields, oldAttributes, attributes, body)

    handlePostResponse(url, objectType, objectID, json, token)

end

# Internal convenience function.  Deletes an object from the repository.
function deleteRepositoryObject(token::AbstractString,
                              objectType::AbstractString,
                              searchKey::Dict{Symbol, Any}
)

    global verbose

    if token == ""
        throw(ArgumentError("`token' cannot be empty"))
    end

    objectID = get(searchKey, :id, 0)
    name = get(searchKey, :name, "")

    object = getObjectInfo(token, objectType, objectID, name)

    # If objectID was not supplied, it will now be available
    objectID = get(object, "id", 0)

    local url = ObjectEndpoint * "/" * objectType * "/$(objectID)"
    local debugID = "(all)"

    if verbose
        println("DELETE $url")
        println("X-Auth-Token: $token")
    end

    resp = Requests.delete(url,
        headers = Dict("X-Auth-Token" => token, "Content-type" => "application/json")
    )

    if statuscode(resp) != 204
        error("Error deleting $(objectType), id = $(objectID).")
    end

    return resp

end

"""
Gets a mapping of custom metric names to Snowflake field names from domain XML.  This list also includes valid dates.

#### Arguments
$(mPulseAPI.readdocs("NodeContent-body"))

#### Returns
$(mPulseAPI.readdocs("CustomMetricMap-structure"))

#### Throws
$(mPulseAPI.readdocs("NodeContent-throws"))
"""
function getCustomMetricMap(body::Any)
    custom_metrics = Dict()

    cmets = getXMLNode(body, "CustomMetrics")

    if !isa(cmets, XMLElement)
        return custom_metrics
    end

    for node in child_elements(cmets)
        attributes = attributes_dict(node)
        if attributes["inactive"] == "false"
            custom_metric = Dict(
                "index" => parse(Int, attributes["index"], 10),
                "fieldname" => "custommetric" * attributes["index"],
                "lastModified" => iso8601ToDateTime(attributes["lastModified"]),
                "description" => attributes["description"]
            )

            datatypeNode = getXMLNode(node, "DataType")
            if isa(datatypeNode, XMLElement)
                custom_metric["dataType"] = attributes_dict(datatypeNode)
            end

            colorNode = getXMLNode(node, "MetricColors")
            if isa(colorNode, XMLElement)
                colors = attributes_dict(colorNode)
                custom_metric["colors"] = map(k->colors[k], sort(collect(keys(colors))))
            end

            custom_metrics[attributes["name"]] = custom_metric
        end
    end

    return custom_metrics
end

"""
Gets a mapping of custom timer names to Snowflake field names from domain XML.  This list also includes valid dates.

#### Arguments
$(mPulseAPI.readdocs("NodeContent-body"))

#### Returns
$(mPulseAPI.readdocs("CustomTimerMap-structure"))

#### Throws
$(mPulseAPI.readdocs("NodeContent-throws"))
"""
function getCustomTimerMap(body::Any)
    custom_timers = Dict()

    ctims = getXMLNode(body, "CustomTimers")

    if !isa(ctims, XMLElement)
        return custom_timers
    end

    for node in child_elements(ctims)
        attributes = attributes_dict(node)
        if attributes["inactive"] == "false"
            custom_timer = Dict(
                "index" => parse(Int, attributes["index"], 10),
                "fieldname" => "customtimer" * attributes["index"],
                "mpulseapiname" => "CustomTimer" * attributes["index"],
                "lastModified" => iso8601ToDateTime(attributes["lastModified"]),
                "description" => attributes["description"]
            )

            colorNodes = getXMLNode(node, "TimingColors")
            if isa(colorNodes, XMLElement)
                colors = Dict[]
                for colorNode in child_elements(colorNodes)
                    color = attributes_dict(colorNode)
                    push!(colors, color)
                end
                custom_timer["colors"] = colors
            end

            custom_timers[attributes["name"]] = custom_timer
        end
    end

    return custom_timers
end

"""
Gets the content of a node

#### Arguments
$(mPulseAPI.readdocs("NodeContent-body"))

`nodeName::AbstractString`
:    The node whose contents shoudl be returned

`default::Any`
:    A default value to return if the required node was not found

#### Returns
`{AbstractString|Number|Boolean}` The content of the requested node cast to the same type as `default` or the value of `default` if the node was not found

#### Throws
$(mPulseAPI.readdocs("NodeContent-throws"))
"""
function getNodeContent(body::Any, nodeName::AbstractString, default::Any)
    node = getXMLNode(body, nodeName)

    # If we have a valid node, get its textContent
    if isa(node, LightXML.XMLElement)
        value = content(node)

        # If the default value passed in was a Number, then we cast value to a Number
        if isa(default, AbstractFloat)
            value = parse(Float64, value, 10)
        elseif isa(default, Int)
            value = parse(Int, value, 10)
        elseif isa(default, Bool)
            if lowercase(value) == "true"
                value = true
            else
                value = false
            end
        end
    else
        value = default
    end

    return value
end

# Internal convenience function
function getXMLNode(body::Any, nodeName::AbstractString)
    if isa(body, AbstractString)
        xdoc = parse_string(body)
        xroot = root(xdoc)
    elseif isa(body, LightXML.XMLElement)
        xroot = body
    elseif isa(body, Dict) && haskey(body, "body")
        xroot = body["body"]
    else
        throw(ArgumentError("bodyXML must either be an XML String, a LightXML.XMLElement or a Dict() with a `body` element. $(typeof(body)) is unknown."))
    end

    return find_element(xroot, nodeName)
end

# Internal convenience function used to POST to repository
function postHttpRequest(url::AbstractString, objectType::AbstractString, objectID::Int64, json::Dict{AbstractString, Any}, token::AbstractString)
    resp = nothing

    try
        resp = Requests.post(url,
                             json = json,
                             headers = Dict("X-Auth-Token" => token, "Content-type" => "application/json")
                        )
    catch er
        if isa(er, Base.UVError)
            error("TCP timeout")
        else
            error("We have not encountered this error before.  Please report this. Timestamp: $(round(Int, datetime2unix(now())))")
        end
    end

    # 400 - Bad request.  The URL or JSON is invalid
    # 404 - Not found.  The requested object does not exist
    if statuscode(resp) == 400 || statuscode(resp) == 404
        throw(mPulseAPIException("Error updating $(objectType) $(objectID)", resp))
    end

    return resp
end

# Internal convenience function for handling POST REST API Responses
function handlePostResponse(url::AbstractString, objectType::AbstractString, objectID::Int64, json::Dict{AbstractString, Any}, token::AbstractString)
    count = 0
    while count <= 5
        resp = postHttpRequest(url, objectType, objectID, json, token)

        if statuscode(resp) == 204 # Success
            return resp
        elseif statuscode(resp) == 401 # Unauthorized.  The security token is missing or invalid.
            # Retry once
            if count > 1
                throw(mPulseAPIAuthException(resp))
            end

            count += 1
        else # Internal server error.  Try again later. Expecting 500 < resp < 509
            # Retry up to 5 times
            if count <= 5
                count += 1
            else
                throw(mPulseAPIBugException(resp))
            end
        end
    end
end

# Internal convenience function for building object JSON entry used in POST
function buildPostJSON(
    objectType::AbstractString,
    objectID::Int64,
    objectFields::Dict=Dict(),
    oldAttributes::Dict=Dict(),
    attributes::Dict=Dict(),
    body::Union{AbstractString, LightXML.XMLElement}=""
)
    # Initialize JSON Dict
    json = Dict{AbstractString, Any}()
    json["type"] = objectType
    json["id"] = objectID

    # If attributes is supplied, update the object’s attributes field
    if !isempty(attributes)
        attributes = convert(Dict{AbstractString, Any}, attributes)

        if objectType == "statisticalmodel"
        # Merge existing and new attributes needed for statisticalmodel
            for key in keys(oldAttributes)
                if !haskey(attributes, key)
                    attributes[key] = oldAttributes[key]
                end
            end
        end

        attributesDict = []

        for (key, val) in attributes
            push!(attributesDict, Dict("name" => key, "value"=> val))
        end

        json["attributes"] = attributesDict

    end

    # If any objectFields are supplied by the user, update these in the object (if it exists)
    if !isempty(objectFields)
        for (key, val) in objectFields
            json[key] = val
        end
    end

    # If the bodyError argument is supplied, update this in the object
    if body != ""
        if isa(body, AbstractString)
            try
                xdoc = parse_string(body)
                xroot = root(xdoc)
            catch
                if objectType == "alert" || objectType == "statisticalmodel"
                    error("errorXML is not formatted correctly")
                else
                    error("body keyword argument is not formatted correctly")
                end
            end
        else
            body = string(body)
        end
        json["body"] = body
    end

    return json
end

# Internal convenience function used to GET from repository
function getHttpRequest(token::AbstractString, objectType::AbstractString, searchKey::Dict{Symbol, Any}, isKeySet::Bool)

    local url = ObjectEndpoint * "/" * objectType
    local query = Dict()
    local debugID = "(all)"

    # Adjust query URL to use ID or search attribute depending on which is passed in
    for (k, v) in searchKey
        if isa(v, Number) ? v > 0 : v != ""
            if k == :id
                url *= "/$(v)"
            else
                query[k] = v
            end
            debugID = "$(k)=$(v)"
            break
        end
    end

    if verbose
        println("GET $url")
        println("X-Auth-Token: $token")
        println(query)
    end

    resp = Requests.get(url, headers=Dict("X-Auth-Token" => token), query=query)

    respStatusCode = statuscode(resp)

    if respStatusCode == 401
        throw(mPulseAPIAuthException(resp))
    elseif respStatusCode == 500
        throw(mPulseAPIBugException(resp))
    elseif respStatusCode != 200
        throw(mPulseAPIException("Error fetching $(objectType) $(debugID)", resp))
    end

    # Do not use Requests.json as that expects UTF-8 data, and mPulse API's response is ISO-8859-1
    json = join(map(Char, resp.data))
    object = JSON.parse(json)

      # If calling by a searchKey other than ID, the return value will be a Dict with a single key="objects"
    # and value set to an array of domain objects.
    # If searching by ID, then the object is returned, so turn it into an array
    if haskey(object, "objects") && isa(object["objects"], Array)
        object_list = object["objects"]
    else
        object_list = [object]
    end

    if length(object_list) == 0
        throw(mPulseAPIException("An object matching $debugID was not returned", resp))
    # If caller has passed in a filter key, then we should only get a single object
    elseif isKeySet && length(object_list) > 1
        throw(mPulseAPIException("Found too many matching objects with IDs=($(map(d->d["id"], object_list)))", resp))
    end

    for object in object_list
        if !isempty(object["body"])
            # Convert body string to an actual XML root object
            xdoc = parse_string(object["body"])
            xroot = root(xdoc)
        else
            # LightXML cannot handle empty strings, so just fake it
            xdoc = XMLDocument()
            xroot = create_root(xdoc, "")
        end

        object["body"] = xroot
        object["created"] = iso8601ToDateTime(object["created"])
        object["lastModified"] = iso8601ToDateTime(object["lastModified"])

        object["attributes"] = Dict(
            map(attr -> attr["name"] => fixJSONDataType(attr["value"]), object["attributes"])
        )

        delete!(object, "schemaVersion")
        delete!(object, "type")
        delete!(object, "effectivePermissions")

        writeObjectToCache(objectType, searchKey, object)

    end

    return object_list

end


text:
#Generic Branch Models
abstract type AbstractBranchFormulation <: AbstractDeviceFormulation end
abstract type AbstractBoundedBranchFormulation <: AbstractBranchFormulation end

#Abstract Line Models
struct StaticLine <: AbstractBranchFormulation end
struct StaticLineBounds <: AbstractBoundedBranchFormulation end
struct StaticLineUnbounded <: AbstractBranchFormulation end
struct FlowMonitoredLine <: AbstractBranchFormulation end

#Abstract Transformer Models
struct StaticTransformer <: AbstractBranchFormulation end
struct StaticTransformerBounds <: AbstractBoundedBranchFormulation end
struct StaticTransformerUnbounded <: AbstractBranchFormulation end

# Not implemented yet
# struct TapControl <: AbstractBranchFormulation end
# struct PhaseControl <: AbstractBranchFormulation end

#################################### Branch Variables ##################################################
# Because of the way we integrate with PowerModels, most of the time PowerSimulations will create variables
# for the branch flows either in AC or DC.
flow_variables!(
    psi_container::PSIContainer,
    ::Type{<:PM.AbstractPowerModel},
    ::IS.FlattenIteratorWrapper{<:PSY.ACBranch},
) = nothing

function flow_variables!(
    psi_container::PSIContainer,
    ::Type{<:StandardPTDFModel},
    devices::IS.FlattenIteratorWrapper{B},
) where {B <: PSY.ACBranch}
    add_variable!(psi_container, devices, make_variable_name(FLOW_ACTIVE_POWER, B), false)
    return
end

#################################### Flow Variable Bounds ##################################################
function _get_constraint_data(
    devices::IS.FlattenIteratorWrapper{B},
) where {B <: PSY.ACBranch}
    constraint_infos = Vector{DeviceRangeConstraintInfo}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        limit_values = (min = -1 * PSY.get_rate(d), max = PSY.get_rate(d))
        name = PSY.get_name(d)
        services_ub = Vector{Symbol}()
        for service in PSY.get_services(d)
            SR = typeof(service)
            push!(services_ub, Symbol("R$(PSY.get_name(service))_$SR"))
        end
        constraint_infos[ix] =
            DeviceRangeConstraintInfo(name, limit_values, services_ub, Vector{Symbol}())
    end
    return constraint_infos
end

function branch_rate_bounds!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{B},
    model::DeviceModel{B, <:AbstractBranchFormulation},
    ::Type{<:PM.AbstractDCPModel},
) where {B <: PSY.ACBranch}
    constraint_infos = _get_constraint_data(devices)
    set_variable_bounds!(psi_container, constraint_infos, FLOW_ACTIVE_POWER, B)
    return
end

function branch_rate_bounds!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{B},
    model::DeviceModel{B, <:AbstractBranchFormulation},
    ::Type{<:PM.AbstractPowerModel},
) where {B <: PSY.ACBranch}
    constraint_infos = _get_constraint_data(devices)
    set_variable_bounds!(psi_container, constraint_infos, FLOW_ACTIVE_POWER_FROM_TO, B)
    set_variable_bounds!(psi_container, constraint_infos, FLOW_ACTIVE_POWER_TO_FROM, B)
    return
end

#################################### Rate Limits constraint_infos ###############################
function branch_rate_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{B},
    model::DeviceModel{B, <:AbstractBranchFormulation},
    ::Type{<:PM.AbstractActivePowerModel},
    feedforward::Nothing,
) where {B <: PSY.ACBranch}
    constraint_infos = _get_constraint_data(devices)
    device_range!(
        psi_container,
        RangeConstraintSpecInternal(
            constraint_infos,
            make_constraint_name(RATE_LIMIT, B),
            make_variable_name(FLOW_ACTIVE_POWER, B),
        ),
    )
    return
end

function branch_rate_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{B},
    model::DeviceModel{B, <:AbstractBranchFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Nothing,
) where {B <: PSY.ACBranch}
    range_data = [(PSY.get_name(h), PSY.get_rate(h)) for h in devices]
    rating_constraint!(
        psi_container,
        range_data,
        make_constraint_name(RATE_LIMIT_FT, B),
        (
            make_variable_name(FLOW_ACTIVE_POWER_FROM_TO, B),
            make_variable_name(FLOW_REACTIVE_POWER_FROM_TO, B),
        ),
    )

    rating_constraint!(
        psi_container,
        range_data,
        make_constraint_name(RATE_LIMIT_TF, B),
        (
            make_variable_name(FLOW_ACTIVE_POWER_TO_FROM, B),
            make_variable_name(FLOW_REACTIVE_POWER_TO_FROM, B),
        ),
    )
    return
end

############################## Flow Limits Constraints #####################################
function branch_flow_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{PSY.MonitoredLine},
    model::DeviceModel{PSY.MonitoredLine, FlowMonitoredLine},
    ::Type{T},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {T <: PM.DCPPowerModel}
    constraint_infos = Vector{PSI.DeviceRangeConstraintInfo}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        if PSY.get_flow_limits(d).to_from != PSY.get_flow_limits(d).from_to
            @info("Flow limits in Line $(PSY.get_name(d)) aren't equal. The minimum will be used in formulation $(T)")
        end
        limit = min(
            PSY.get_rate(d),
            PSY.get_flow_limits(d).to_from,
            PSY.get_flow_limits(d).from_to,
        )
        minmax = (min = -1 * limit, max = limit)
        constraint_infos[ix] = DeviceRangeConstraintInfo(PSY.get_name(d), minmax)
    end
    device_range!(
        psi_container,
        RangeConstraintSpecInternal(
            constraint_infos,
            make_constraint_name(FLOW_LIMIT, PSY.MonitoredLine),
            make_variable_name(FLOW_ACTIVE_POWER, PSY.MonitoredLine),
        ),
    )
    return
end

function branch_flow_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{PSY.MonitoredLine},
    model::DeviceModel{PSY.MonitoredLine, FlowMonitoredLine},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
)
    names = Vector{String}(undef, length(devices))
    limit_values_FT = Vector{MinMax}(undef, length(devices))
    limit_values_TF = Vector{MinMax}(undef, length(devices))
    to = Vector{PSI.DeviceRangeConstraintInfo}(undef, length(devices))
    from = Vector{PSI.DeviceRangeConstraintInfo}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        limit_values_FT[ix] = (
            min = -1 * PSY.get_flow_limits(d).from_to,
            max = PSY.get_flow_limits(d).from_to,
        )
        limit_values_TF[ix] = (
            min = -1 * PSY.get_flow_limits(d).to_from,
            max = PSY.get_flow_limits(d).to_from,
        )
        names[ix] = PSY.get_name(d)
        to[ix] = DeviceRangeConstraintInfo(names[ix], limit_values_FT[ix])
        from[ix] = DeviceRangeConstraintInfo(names[ix], limit_values_TF[ix])
    end

    device_range!(
        psi_container,
        RangeConstraintSpecInternal(
            to,
            make_constraint_name(FLOW_LIMIT_FROM_TO, PSY.MonitoredLine),
            make_variable_name(FLOW_ACTIVE_POWER_FROM_TO, PSY.MonitoredLine),
        ),
    )
    device_range!(
        psi_container,
        RangeConstraintSpecInternal(
            from,
            make_constraint_name(FLOW_LIMIT_TO_FROM, PSY.MonitoredLine),
            make_variable_name(FLOW_ACTIVE_POWER_TO_FROM, PSY.MonitoredLine),
        ),
    )
    return
end


text:
#Tools to help get a feel for Tr0
module Tr0Tools

#==Aliases
===============================================================================#
const DataWord = UInt32

#==Constants
===============================================================================#
const WRITEBLOCK_SYNCWORD = DataWord(0x4)

#==Main Types
===============================================================================#
struct BlockHeader
	_type::DataWord
	_size::DataWord
end

#==Helper Functions
===============================================================================#
function corruptword_exception(io::IO, w::DataWord, expected::DataWord)
	pos = position(io) - sizeof(DataWord)
	pos = hex(pos)
	w = hex(w)
	expected = hex(expected)
	return "Corrupt word @ 0x$pos: 0x$w, 0x$expected"
end

function readsyncword(io::IO)
	w = read(io, DataWord)
	if w != WRITEBLOCK_SYNCWORD
		throw(corruptword_exception(io, w, WRITEBLOCK_SYNCWORD))
	end
end

#Data read:
function _dread(io::IO, ::Type{BlockHeader})
	readsyncword(io)
	_type = read(io, DataWord)
	readsyncword(io)
	_size = read(io, DataWord)
	return BlockHeader(_type, _size)
end

function _show(io::IO, hdr::BlockHeader, pos::Int)
	print(io, "Block: 0x", hex(WRITEBLOCK_SYNCWORD))
	print(io, " 0x", hex(hdr._type))
	print(io, " 0x", hex(WRITEBLOCK_SYNCWORD))
	print(io, " 0x", hex(hdr._size))
	print(io, " (start 0x", hex(pos), ")")
	println(io)
end

#==Main functions
===============================================================================#
function dumpsegments(io::IO, filepath::String)
	r = open(filepath)
	blockcount = 0
	totalsize = 0

	while !eof(r)
		pos = position(r)
		hdr = _dread(r, BlockHeader)
			_show(io, hdr, pos)
			blockcount += 1
			totalsize += hdr._size

		seek(r, position(r)+hdr._size)
		blksize = read(r, DataWord)
		if blksize != hdr._size
			throw(corruptword_exception(r, blksize, hdr._size))
		end
	end

	close(r)
	println(io, "Blocks read: $blockcount, total size: $totalsize.")
end

dumpsegments(filepath::String) = dumpsegments(STDOUT, filepath)

end #module


text:
#===============================================================================

    CalibrationParameters.jl

    macro-generated struct of parameters
    for use in calibration analysis

    Author: Tom Price
    Date:   July 2019

===============================================================================#

#===============================================================================
    defaults >>
===============================================================================#

const DEFAULT_CAL_DCV               = true
const DEFAULT_CAL_DYE_IN            = :FAM
const DEFAULT_CAL_DYES_TO_FILL      = Vector{Symbol}()
const DEFAULT_NORM_SUBTRACT_WATER   = false
const DEFAULT_DCV_K_METHOD          = well_proc_vec

#===============================================================================
    field definitions >>
===============================================================================#

const CAL_ARGS_FIELD_DEFS = [
    Field(:dcv,                     Bool,               DEFAULT_CAL_DCV),
    Field(:dye_in,                  Symbol,             Meta.quot(DEFAULT_CAL_DYE_IN)),
    Field(:dyes_to_fill,            Vector{Symbol},     DEFAULT_CAL_DYES_TO_FILL),
    Field(:subtract_water,          Bool,               DEFAULT_NORM_SUBTRACT_WATER),
    Field(:k_method,                KMethod,            DEFAULT_DCV_K_METHOD)]

#===============================================================================
    macro calls >>
===============================================================================#

## generate struct and constructor
SCHEMA = CAL_ARGS_FIELD_DEFS
# println(@macroexpand @make_struct_from_SCHEMA CalibrationParameters)
# println(@macroexpand @make_constructor_from_SCHEMA CalibrationParameters)
@make_struct_from_SCHEMA CalibrationParameters
@make_constructor_from_SCHEMA CalibrationParameters

#===============================================================================
    constant >>
===============================================================================#

const DEFAULT_CAL_ARGS = CalibrationParameters()

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
using Test
using Flux: onehotbatch, mse, crossentropy, logitcrossentropy,
            σ, binarycrossentropy, logitbinarycrossentropy

const ϵ = 1e-7

@testset "losses" begin
  # First, regression-style y's
  y = [1, 1, 0, 0]
  ŷ = [.9, .1, .1, .9]

  @testset "mse" begin
    @test mse(ŷ, y) ≈ (.1^2 + .9^2)/2
  end

  # Now onehot y's
  y = onehotbatch([1, 1, 0, 0], 0:1)
  ŷ = [.1 .9; .9 .1; .9 .1; .1 .9]'
  v = log(.1 / .9)
  logŷ = [v 0.0; 0.0 v; 0.0 v; v 0.0]'
  lossvalue = 1.203972804325936

  @testset "crossentropy" begin
    @test crossentropy(ŷ, y) ≈ lossvalue
  end

  @testset "logitcrossentropy" begin
    @test logitcrossentropy(logŷ, y) ≈ lossvalue
  end

  @testset "weighted_crossentropy" begin
    @test crossentropy(ŷ, y, weight = ones(2)) ≈ lossvalue
    @test crossentropy(ŷ, y, weight = [.5, .5]) ≈ lossvalue/2
    @test crossentropy(ŷ, y, weight = [2, .5]) ≈ 1.5049660054074199
  end

  @testset "weighted_logitcrossentropy" begin
    @test logitcrossentropy(logŷ, y, weight = ones(2)) ≈ lossvalue
    @test logitcrossentropy(logŷ, y, weight = [.5, .5]) ≈ lossvalue/2
    @test logitcrossentropy(logŷ, y, weight = [2, .5]) ≈ 1.5049660054074199
  end

  logŷ, y = randn(3), rand(3)
  @testset "binarycrossentropy" begin
    @test binarycrossentropy.(σ.(logŷ), y; ϵ=0) ≈ -y.*log.(σ.(logŷ)) - (1 .- y).*log.(1 .- σ.(logŷ))
    @test binarycrossentropy.(σ.(logŷ), y) ≈ -y.*log.(σ.(logŷ) .+ eps.(σ.(logŷ))) - (1 .- y).*log.(1 .- σ.(logŷ) .+ eps.(σ.(logŷ)))
  end

  @testset "logitbinarycrossentropy" begin
    @test logitbinarycrossentropy.(logŷ, y) ≈ binarycrossentropy.(σ.(logŷ), y; ϵ=0)
  end
end


text:
using Documenter, BioSimulator

# setup GR backend for Travis CI
# ENV["GKSwstype"] = "100"
# ENV["PLOTS_TEST"] = "true"

makedocs(
  doctest  = false,
  format   = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
  modules  = [BioSimulator],
  clean    = true,
  sitename = "BioSimulator.jl",
  authors  = "Alfonso Landeros, Timothy Stutz, Kevin L. Keys, Mary E. Sehl, Alexander Alekseyenko, Janet S. Sinsheimer, Kenneth Lange",
  pages = [
    "Home" => "index.md",
    "Model Specification" => "man/model-specification.md",
    "Simulation" => "man/simulation.md",
    "Plotting" => "man/plotting.md",
    "Tutorials" => "man/tutorials.md",
    "Algorithms and References" => "man/algorithms.md",
    "Citing" => "man/citing.md"
  ]
)

deploydocs(
  repo   = "github.com/alanderos91/BioSimulator.jl.git",
)


text:
# test two properties of oscillating matrix
# For an n x n oscillating matrix A, the following holds:
# 1. A has n distinct and postive eigenvalues λ_1 > λ_2 > ... > λ_n > 0
# 2. The ith eigenvector, corresponding to λ_i in the above ordering has
#    exactly i-1 sign changes.

mode = rand(1:2) # 2 different modes
n = rand(2:10)
A = matrixdepot("oscillate", n, mode)
try
    matrixdepot("oscillate", n, 4)
catch ArgumentError
    println("oscillate invalid mode value")
end

eva, evc = eig(A)
for i in eva
    @test i > 0
end

p = sortperm(eva, rev = true)
evc = evc[:,p]

# compute the number of sign change
function num_sign_change(v)
    signv = sign(v)
    change = signv[1]
    num = 0
    for i in signv
        if i != change
            num += 1
            change = i
        end
    end
    return num
end

k = rand(1:n)
@test num_sign_change(evc[:,k]) == k - 1

A = matrixdepot("oscillate", n)

println("'oscillate' passed test...")


text:
# ============================================================================ #
# DictSet
# ============================================================================ #
# TODO: use IntSet? - benchmark!
"""
Basic dict-like thing that holds multiple values for each key in a set.
"""
struct DictSet{K}
    data::Dict{K, IntSet}
end

DictSet{K}() where K = DictSet{K}(Dict{K, IntSet}())

function Base.getindex(ds::DictSet{K}, k::K) where K
    if haskey(ds.data, k)
        ds.data[k]
    else
        IntSet()
    end
end

function Base.insert!(ds::DictSet{K}, k::K, v::Int) where K
    if haskey(ds.data, k)
        push!(ds.data[k], v)
    else
        ds.data[k] = IntSet(v)
    end
    v
end

function Base.delete!(ds::DictSet{K}, k::K) where K
    delete!(ds.data, k)
    ds
end

function Base.delete!(ds::DictSet{K}, k::K, v::Int) where K
    if haskey(ds.data, k)
        s = ds.data[k]
        if length(s) > 1
            delete!(s, v)
        else
            delete!(ds.data, k)
        end
    end
    ds
end

function Base.delete!(ds::DictSet{K}, vals::IntSet) where K
    for (k, v) in ds.data
        setdiff!(v, vals)
        if isempty(v)
            delete!(ds.data, k)
        end
    end
end

Base.keys(ds::DictSet) = keys(ds.data)
Base.values(ds::DictSet) = union(values(ds.data)...)
Base.haskey(ds::DictSet) = haskey(ds.data)
Base.isempty(ds::DictSet) = isempty(ds.data)

# ============================================================================ #
# THE SAL TYPE
# ============================================================================ #
mutable struct SAL{T}
    nodelist::Vector{DictSet{Int}}

    # TODO: store data
    # reverse dict for maximal simplices?
    lbl_vals::IntSet
    lbl_holes::IntSet
end

SAL{T}() where {T} = SAL{T}(DictSet{Int}[], IntSet(), IntSet())
SAL() = SAL{Float64}()

Base.show(io::IO, sal::SAL) = print(io, "SAL with $(length(sal)) nodes")

# Pretty printing.
function Base.print(io::IO, sal::SAL)
    println("SAL:")
    for i in eachindex(sal)
        print("$(lpad(i, 5)): [")
        for k in sort(collect(keys(sal.nodelist[i].data)))
            for v in sort(collect(sal.nodelist[i].data[k]))
                print("($i, $(iszero(k) ? "φ" : k), $v)")
            end
        end
        println("]")
    end
end

# ============================================================================ #
# GETTERS, SETTERS, UTILS
# ============================================================================ #
function getlabel!(sal::SAL)
    if isempty(sal.lbl_holes)
        lbl = length(sal.lbl_vals) + 1
    else
        lbl = first(sal.lbl_holes)
        delete!(sal.lbl_holes, lbl)
    end
    push!(sal.lbl_vals, lbl)
    lbl
end

Base.length(sal::SAL) = length(sal.nodelist)
Base.eachindex(sal::SAL) = 1:length(sal.nodelist)
Base.isempty(sal::SAL) = all(isempty.(sal.nodelist))

function insnode!(sal::SAL, cur::Int, nxt::Int, lbl::Int)
    while length(sal) < cur
        push!(sal.nodelist, DictSet{Int}())
    end

    insert!(sal.nodelist[cur], nxt, lbl)
end

function labels(sal::SAL, from::Int, to::Int)
    sal.nodelist[from][to]
end

function alllabels(sal::SAL)
    sal.lbl_values
end

function maximalsimplices(sal::SAL)
    maxsxs = Dict{Int, Vector{Int}}()
    for l in sal.lbl_vals
        maxsxs[l] = []
    end
    for i in eachindex(sal)
        for v in values(sal.nodelist[i])
            push!(maxsxs[v], i)
        end
    end
    maxsxs
end

"""
Check if σ1 is a face of σ2. Assumes the simplices are properly formatted.
"""
function isface(σ1, σ2)
    i = 1; j = 1
    while i <= length(σ1) && j <= length(σ2)
        if σ1[i] > σ2[j]
            j += 1
        elseif σ1[i] == σ2[j]
            i += 1; j += 1
        else
            return false
        end
    end
    i == length(σ1) + 1
end

# ============================================================================ #
# OPERATIONS
# ============================================================================ #
# A σ is in the complex if it is contained in one of the maximal sx's.
function Base.in(σ, sal::SAL)
    σ = check_simplex(σ)
    !isempty(inmaximals(sal, σ))
end

"""
List maximal simplices that contain σ.
"""
function inmaximals(sal::SAL, σ)
    last(σ) > length(sal) && return IntSet()

    s = copy(labels(sal, σ[1], σ[2]))
    for l in zip(σ[2:end-1], σ[3:end])
        intersect!(s, labels(sal, l...))
        isempty(s) && return s
    end
    s
end

function Base.push!(sal::SAL, σ)
    # TODO: remove faces of simplex!
    # Find maximal simplices that are contained in σ and remove them
    σ = check_simplex(σ)
    σ in sal && return σ

    lbl = getlabel!(sal)

    for (i, li) in enumerate(σ), (j, lj) in enumerate(σ[i+1:end])
        insnode!(sal, li, lj, lbl)
    end
    insnode!(sal, σ[end], 0, lbl)
    σ
end

function deletemaximals!(sal::SAL, lbls::IntSet)
    for i in eachindex(sal)
        delete!(sal.nodelist[i], lbls)
    end
    union!(sal.lbl_holes, lbls)
    setdiff!(sal.lbl_vals, lbls)
end

"""

"""
function facetswithout(σ, τ)
end

function Base.delete!(sal::SAL, σ)
    # TODO: multiple at once!
    # insert facets that do not contain σ
    deleted_lbls = IntSet()
    deleted_sxs  = Vector{Int}[]
    for (lbl, τ) in maximalsimplices(sal)
        if isface(τ, σ)
            push!(delted_lbls, lbl)
            push!(delted_sxs, τ)
        end
    end
    deletemaximals!(sal, deleted_lbls)
    for τ in delted_sxs
        for υ in facetswithout(τ, σ)
            insert!(sal, υ)
        end
    end

    # vsi k vsebujejo σ = greš po drevesu kot σ in delaš intersectione
    for (i, li) in enumerate(σ)
    end
end


text:
# Initialization algorithms
#
#   Each algorithm is represented by a subtype of SeedingAlgorithm
#
#   Let alg be an instance of such an algorithm, then it should
#   support the following usage:
#
#       initseeds!(iseeds, alg, X)
#       initseeds_by_costs!(iseeds, alg, costs)
#
#   Here:
#       - iseeds:   a vector of resultant indexes of the chosen seeds
#       - alg:      the seeding algorithm instance
#       - X:        the data matrix, each column being a sample
#       - costs:    pre-computed pairwise cost matrix.
#
#   This function returns iseeds
#

@compat abstract type SeedingAlgorithm end

initseeds(alg::SeedingAlgorithm, X::RealMatrix, k::Integer) =
    initseeds!(Vector{Int}(k), alg, X)

initseeds_by_costs(alg::SeedingAlgorithm, costs::RealMatrix, k::Integer) =
    initseeds_by_costs!(Vector{Int}(k), alg, costs)

seeding_algorithm(s::Symbol) =
    s == :rand ? RandSeedAlg() :
    s == :kmpp ? KmppAlg() :
    s == :kmcen ? KmCentralityAlg() :
    error("Unknown seeding algorithm $s")

initseeds(algname::Symbol, X::RealMatrix, k::Integer) =
    initseeds(seeding_algorithm(algname), X, k)::Vector{Int}

initseeds_by_costs(algname::Symbol, costs::RealMatrix, k::Integer) =
    initseeds_by_costs(seeding_algorithm(algname), costs, k)

initseeds(iseeds::Vector{Int}, X::RealMatrix, k::Integer) = iseeds
initseeds_by_costs(iseeds::Vector{Int}, costs::RealMatrix, k::Integer) = iseeds

function copyseeds!(S::DenseMatrix, X::DenseMatrix, iseeds::AbstractVector)
    d = size(X, 1)
    n = size(X, 2)
    k = length(iseeds)
    (size(X,1) == d && size(S) == (d, k)) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))

    for j = 1:k
        copy!(view(S,:,j), view(X,:,iseeds[j]))
    end
    return S
end

copyseeds{T}(X::DenseMatrix{T}, iseeds::AbstractVector) =
    copyseeds!(Matrix{T}(size(X,1), length(iseeds)), X, iseeds)

function check_seeding_args(n::Integer, k::Integer)
    k >= 1 || error("The number of seeds must be positive.")
    k <= n || error("Attempted to select more seeds than samples.")
end

# Random seeding
#
#   choose an arbitrary subset as seeds
#

type RandSeedAlg <: SeedingAlgorithm end

initseeds!(iseeds::IntegerVector, alg::RandSeedAlg, X::RealMatrix) =
    sample!(1:size(X,2), iseeds; replace=false)

initseeds_by_costs!(iseeds::IntegerVector, alg::RandSeedAlg, X::RealMatrix) =
    sample!(1:size(X,2), iseeds; replace=false)

# Kmeans++ seeding
#
#   D. Arthur and S. Vassilvitskii (2007).
#   k-means++: the advantages of careful seeding.
#   18th Annual ACM-SIAM symposium on Discrete algorithms, 2007.
#

type KmppAlg <: SeedingAlgorithm end

function initseeds!(iseeds::IntegerVector, alg::KmppAlg, X::RealMatrix, metric::PreMetric)
    n = size(X, 2)
    k = length(iseeds)
    check_seeding_args(n, k)

    # randomly pick the first center
    p = rand(1:n)
    iseeds[1] = p

    if k > 1
        mincosts = colwise(metric, X, view(X,:,p))
        mincosts[p] = 0

        # pick remaining (with a chance proportional to mincosts)
        tmpcosts = zeros(n)
        for j = 2:k
            p = wsample(1:n, mincosts)
            iseeds[j] = p

            # update mincosts
            c = view(X,:,p)
            colwise!(tmpcosts, metric, X, view(X,:,p))
            updatemin!(mincosts, tmpcosts)
            mincosts[p] = 0
        end
    end

    return iseeds
end

initseeds!(iseeds::IntegerVector, alg::KmppAlg, X::RealMatrix) =
    initseeds!(iseeds, alg, X, SqEuclidean())

function initseeds_by_costs!(iseeds::IntegerVector, alg::KmppAlg, costs::RealMatrix)
    n = size(costs, 1)
    k = length(iseeds)
    check_seeding_args(n, k)

    # randomly pick the first center
    p = rand(1:n)
    iseeds[1] = p

    if k > 1
        mincosts = copy(view(costs,:,p))
        mincosts[p] = 0

        # pick remaining (with a chance proportional to mincosts)
        for j = 2:k
            p = wsample(1:n, mincosts)
            iseeds[j] = p

            # update mincosts
            updatemin!(mincosts, view(costs,:,p))
            mincosts[p] = 0
        end
    end

    return iseeds
end

kmpp(X::RealMatrix, k::Int) = initseeds(KmppAlg(), X, k)
kmpp_by_costs(costs::RealMatrix, k::Int) = initseeds(KmppAlg(), costs, k)

# K-medoids initialization based on centrality
#
#   Hae-Sang Park and Chi-Hyuck Jun.
#   A simple and fast algorithm for K-medoids clustering.
#   doi:10.1016/j.eswa.2008.01.039
#

type KmCentralityAlg <: SeedingAlgorithm end

function initseeds_by_costs!(iseeds::IntegerVector, alg::KmCentralityAlg, costs::RealMatrix)
    n = size(costs, 1)
    k = length(iseeds)
    k <= n || error("Attempted to select more seeds than samples.")

    # compute score for each item
    coefs = vec(sum(costs, 2))
    for i = 1:n
        @inbounds coefs[i] = inv(coefs[i])
    end

    # scores[j] = \sum_j costs[i,j] / (\sum_{j'} costs[i,j'])
    #           = costs[i,j] * coefs[i]
    #
    # So this is matrix-vector multiplication
    scores = costs'coefs

    # lower score indicates better seeds
    sp = sortperm(scores)
    for i = 1:k
        @inbounds iseeds[i] = sp[i]
    end
    return iseeds
end

initseeds!(iseeds::IntegerVector, alg::KmCentralityAlg, X::RealMatrix, metric::PreMetric) =
    initseeds_by_costs!(iseeds, alg, pairwise(metric, X))

initseeds!(iseeds::IntegerVector, alg::KmCentralityAlg, X::RealMatrix) =
    initseeds!(iseeds, alg, X, SqEuclidean())


text:
@testset "Row operations" begin
    t = Table(a = [1,2,3])

    @test vcat(t,t)::Table == Table(a = [1,2,3,1,2,3])
    @test empty(t)::Table == Table(a = Int[])

    @test (empty!(t); t == Table(a = Int[]))
    @test (push!(t, (a=2,)); t == Table(a = [2]))
    @test (pushfirst!(t, (a=1,)); t == Table(a = [1,2]))
    @test (append!(t, Table(a=[3])); t == Table(a = [1,2,3]))
    @test (prepend!(t, copy(t)); t == Table(a = [1,2,3,1,2,3]))
    @test (deleteat!(t, 2); t == Table(a = [1,3,1,2,3]))
    @test (insert!(t, 2, (a=0,)); t == Table(a = [1,0,3,1,2,3]))

    @test splice!(t, 2, (a=-1,)) == (a=0,)
    @test t == Table(a = [1,-1,3,1,2,3])
    @test splice!(t, 2, Table(a = [-2])) == (a=-1,)
    @test t == Table(a = [1,-2,3,1,2,3])
    @test splice!(t, 2) == (a=-2,)
    @test t == Table(a = [1,3,1,2,3])

    @test splice!(t, 2:1, (a=-1,)) == Table(a = Int[])
    @test t == Table(a = [1,-1,3,1,2,3])
    @test splice!(t, 2:2, Table(a = [-2])) == Table(a = [-1])
    @test t == Table(a = [1,-2,3,1,2,3])
    @test splice!(t, 1:3) == Table(a=[1,-2,3])
    @test t == Table(a = [1,2,3])

    @test (resize!(t, 10); length(t) == 10)

end


text:
#
# abstract interfaces
#

abstract type AbstractCurve end

# Required methods
"""
	point(C::AbstractCurve,t::Real)
Find the point on curve `C` at parameter value `t`, which should lie in the interval [0,1].
"""
point(c::AbstractCurve,t::Real) = @error "No point() method defined for type $(typeof(c))"

"""
	tangent(C::AbstractCurve,t::Real)
Find the complex number representing the tangent to curve `C` at parameter value `t` in [0,1].
"""
tangent(C::AbstractCurve,t::Real) = @error "No tangent() method defined for type $(typeof(C))"

reverse(C::AbstractCurve) = @error "No reverse() method defined for type $(typeof(C))"

"""
	isfinite(C::AbstractCurve)
Return `true` if the curve is bounded in the complex plane (i.e., does not pass through infinity).
"""
isfinite(C::AbstractCurve) = @error "No isfinite() method defined for type $(typeof(C))"

conj(C::AbstractCurve) = @error "No conj() method defined for type $(typeof(C))"
+(C::AbstractCurve,z::Number) = @error "No addition method defined for type $(typeof(C))"
-(C::AbstractCurve) = @error "No negation method defined for type $(typeof(C))"
*(C::AbstractCurve,z::Number) = @error "No multiplication method defined for type $(typeof(C))"
inv(C::AbstractCurve) = @error "No inversion method defined for type $(typeof(C))"

# Default implementations
"""
	point(C::AbstractCurve,t::AbstractArray)

Vectorize the `point` function for curve `C`.
"""
point(c::AbstractCurve,t::AbstractArray{T}) where T<:Real = [point(c,t) for t in t]

"""
	unittangent(C::AbstractCurve,t::Real)
Find the complex number representing the unit tangent to curve `C` at parameter value `t` in [0,1]. For Lines, Segments, and Rays, the `t` argument is optional.
"""
unittangent(C::AbstractCurve,t::Real) = sign(tangent(C,t))

"""
	normal(C::AbstractCurve,t::Real)
Find the unit complex number in the direction of the leftward-pointing normal to curve `C` at parameter value `t` in [0,1].
"""
normal(c::AbstractCurve,t::Real) = 1im*unittangent(c,t)

+(z::Number,C::AbstractCurve) = +(C,z)
-(C::AbstractCurve,z::Number) = C + (-z)
-(z::Number,C::AbstractCurve) = z + (-C)
*(z::Number,C::AbstractCurve) = *(C,z)
/(C::AbstractCurve,z::Number) = C*(1/z)
/(z::Number,C::AbstractCurve) = z*inv(C)

function arclength(C::AbstractCurve,part=[0,1])
	f = t -> abs(tangent(C,t))
	intadapt(f,part...,DEFAULT[:tol])
end

"""
	plotdata(C::AbstractCurve,n=501)

Compute `n` points along the curve `C` suitable to make a plot of it.
"""
plotdata(C::AbstractCurve) = adaptpoints(t->point(C,t),t->unittangent(C,t),0,1)

show(io::IO,C::AbstractCurve) = print(io,"Complex-valued $(typeof(C))")
show(io::IO,::MIME"text/plain",C::AbstractCurve) = print(io,"Complex-valued $(typeof(C))")

# AbstractClosedCurve
abstract type AbstractClosedCurve <: AbstractCurve end

# Default implementations
function winding(C::AbstractClosedCurve,z::Number)
	# Integrate around the curve
	f = t -> imag(tangent(C,t)/(point(C,t)-z))
	w = intadapt(f,0,1,1e-4)
	return round(Int,w/(2π))
end

isinside(z::Number,C::AbstractClosedCurve) = winding(C,z) != 0
isoutside(z::Number,C::AbstractClosedCurve) = winding(C,z) == 0

#
# generic Curve
#

"""
(type) Smooth curve defined by an explicit function of a real paramerter in [0,1].
"""
struct Curve <: AbstractCurve
	point
	tangent
end

"""
	Curve(f)
	Curve(f,a,b)
Construct a `Curve` object from the complex-valued function `f` accepting an argument in the interval [0,1]. If `a` and `b` are given, they are the limits of the parameter in the call to the supplied `f`. However, the resulting object will be defined on [0,1], which is internally scaled to [a,b].

	Curve(f,df[,a,b])
Construct a curve with point location and tangent given by the complex-valued functions `f` and `df`, respectively, optionally with given limits on the parameter.
"""
Curve(f) = Curve(f,t->fdtangent(f,t))
Curve(f,df,a::Real,b::Real) = Curve(t->f(scaleto(a,b,t)),t->df(scaleto(a,b,t)))
Curve(f,a::Real,b::Real) = Curve(t->f(scaleto(a,b,t)))

# Required methods
point(C::Curve,t::Real) = C.point(t)
(C::Curve)(t::Real) = point(C,t)
tangent(C::Curve,t::Real) = C.tangent(t)
conj(C::Curve) = Curve(t->conj(C.point(t)))
reverse(C::Curve) = Curve(t->C.point(1-t))
isfinite(C::Curve) = true

+(C::Curve,z::Number) = Curve(t->C.point(t)+z,C.tangent)
-(C::Curve) = Curve(t->-C.point(t),t->-C.tangent(t))
*(C::Curve,z::Number) = Curve(t->C.point(t)*z,t->C.tangent(t)*z)
inv(C::Curve) = Curve(t->1/C.point(t),t->-C.tangent(t)/C.point(t)^2)

#
# generic ClosedCurve
#

"""
(type) Smooth closed curve defined by an explicit function of a real paramerter in [0,1].
"""
struct ClosedCurve <: AbstractClosedCurve
	curve::Curve
	function ClosedCurve(c::Curve;tol=DEFAULT[:tol])
		@assert isapprox(point(c,0),point(c,1);rtol=tol,atol=tol) "Curve does not close"
		new(c)
	end
end

"""
	ClosedCurve(f; tol=<default>)
	ClosedCurve(f,a,b; tol=<default)
Construct a `ClosedCurve` object from the complex-valued function `point` accepting an argument in the interval [0,1]. The constructor checks whether `f(0)≈f(1)` to tolerance `tol`. If `a` and `b` are given, they are the limits of the parameter in the call to the supplied `f`. However, the resulting object will be defined on [0,1], which is internally scaled to [a,b].

	ClosedCurve(f,df[,a,b]; tol=<default>)
Construct a closed curve with point location and tangent given by the complex-valued functions `f` and `df`, respectively, optionally with given limits on the parameter.
"""
ClosedCurve(f,df=t->fdtangent(f,t);kw...) = ClosedCurve(Curve(f,df;kw...))
ClosedCurve(f,a::Real,b::Real;kw...) = ClosedCurve(Curve(f,a,b;kw...))
ClosedCurve(f,df,a::Real,b::Real;kw...) = ClosedCurve(Curve(f,df,a,b;kw...))

point(C::ClosedCurve,t::Real) = point(C.curve,t)
(C::ClosedCurve)(t::Real) = point(C.curve,t)
tangent(C::ClosedCurve,t::Real) = tangent(C.curve,t)
conj(C::ClosedCurve) = ClosedCurve(conj(C.curve))
reverse(C::ClosedCurve) = ClosedCurve(reverse(C.curve))
isfinite(C::ClosedCurve) = isfinite(C.curve)
+(C::ClosedCurve,z::Number) = ClosedCurve(C.curve+z)
-(C::ClosedCurve) = ClosedCurve(-C.curve)
*(C::ClosedCurve,z::Number) = ClosedCurve(C.curve*z)
inv(C::ClosedCurve) = ClosedCurve(inv(C.curve))

include("lines.jl")
include("rays.jl")
include("segments.jl")
include("circles.jl")
include("arcs.jl")
include("intersections.jl")


text:
"""
Created on Fri 15 May 2020

@author: Yoann Pradat

Performs NMF rank selection by the method of factor stability across sligth perturbations of the V matrix.
"""

using Distances
using Distributed
using Distributions
using DelimitedFiles
using PyCall

@everywhere include("./src/VNMF.jl")
@everywhere using .VNMF

#### folder for saving plots and other files
folder   = "./demo/plot/breast_21"
mkpath(folder)

#### # 1. LOAD DATA
#### # ################################################################################################################

row_types    = vec(readdlm(pwd() * "/../../../data/breast_21/types.txt", ',', String, '\n'))
row_subtypes = vec(readdlm(pwd() * "/../../../data/breast_21/subtypes.txt", ',', String, '\n'))
V = readdlm(pwd() * "/../../../data/breast_21/originalGenomes.txt", ',', Float64, '\n')

# new_order = sortperm(row_types)
# row_types = row_types[new_order]
# row_subtypes = row_subtypes[new_order]
# V = V[new_order, :]
#
# py"""
# import sys
# if "." not in sys.path:
#     sys.path.append(".")
#
# from demo.src.util.plot_VWH import *
#
# def make_plot_V(V, row_types, row_subtypes, filepath):
#     row_groups_colors = {
#         "C>A": "cornflowerblue",
#         "C>G": "lightseagreen",
#         "C>T": "coral",
#         "T>A": "sandybrown",
#         "T>C": "forestgreen",
#         "T>G": "deeppink"
#     }
#
#     fig = plot_V(
#         V                 = V,
#         row_groups        = row_types,
#         row_groups_colors = row_groups_colors,
#         row_labels        = row_subtypes,
#         row_labelsize     = 4
#     )
#     fig.savefig(filepath, bbox_inches="tight")
# """
#
# filepath = joinpath(folder, "plot_V_obs.pdf")
# py"make_plot_V"(V, row_types, row_subtypes, filepath)

#### # 2. SELECT ALGORITHM
#### # ################################################################################################################

#### params for rank selection
rs_params = RSParams(
    K_min  = 1,
    K_max  = 10,
    n_iter = 50,
    pert   = :Multinomial,
    seed   = 123
)

#### nmf algorithm
nmf_global_params = NMFParams(
    init          = :random,
    dist          = Uniform(0, 1),
    max_iter      = 10_000,
    stopping_crit = :conn,
    stopping_iter = 100,
    verbose       = false,
)

nmf_local_params = MUParams(
    β       = 1,
    div     = :β,
    scale_W = false,
    alg     = :mu
)

nmf = NMF(
    solver        = nmf_MU,
    global_params = nmf_global_params,
    local_params  = nmf_local_params
)

#### run rank selection by stability procedure
@time rs_results = rank_by_stability(V, rs_params, nmf)

#### # 3. PLOT RESULTS
#### # ################################################################################################################

#### plot metrics
include(pwd() * "/demo/src/util/plot_metrics.jl")
filepath = joinpath(folder, "rank_stability_metrics.pdf")
plot_metrics_stability(rs_results.df_metrics, filepath)

#### plot factors

py"""
import matplotlib.cm as cm
import sys
if "." not in sys.path:
    sys.path.append(".")

from demo.src.util.plot_VWH import *

def make_plot_W(W_avg, W_std, stab_avg, stab_std, row_types, row_subtypes, filepath):
    row_groups_colors = {
        "C>A": "cornflowerblue",
        "C>G": "lightseagreen",
        "C>T": "coral",
        "T>A": "sandybrown",
        "T>C": "forestgreen",
        "T>G": "deeppink"
    }
    hmap_boundaries = [0, 0.001, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5]
    hmap_cmap = "Blues"

    fig = plot_W_stab(
        W_avg             = W_avg,
        W_std             = W_std,
        stab_avg          = stab_avg,
        stab_std          = stab_std,
        row_groups        = row_types,
        row_groups_colors = row_groups_colors,
        row_labels        = row_subtypes,
        hmap_boundaries   = hmap_boundaries,
        hmap_cmap         = hmap_cmap
    )
    fig.savefig(filepath, bbox_inches="tight")
"""

for K = rs_params.K_min:rs_params.K_max
    filepath = joinpath(folder, "plot_W_rank_$(K).pdf")
    clu_results = rs_results.list_clu_results[K]
    py"make_plot_W"(clu_results.W_avg, clu_results.W_std, clu_results.stab_avg, clu_results.stab_std, row_types, row_subtypes, filepath)
end

#### DEBUG

K = 1

W_init = readdlm(pwd() * "/../../../data/breast_21/W_init_rank$K.txt", ',', Float64, '\n')
H_init = readdlm(pwd() * "/../../../data/breast_21/H_init_rank$K.txt", ',', Float64, '\n')

W_fit = readdlm(pwd() * "/../../../data/breast_21/W_fit_rank$K.txt", ',', Float64, '\n')
H_fit = readdlm(pwd() * "/../../../data/breast_21/H_fit_rank$K.txt", ',', Float64, '\n')

V = readdlm(pwd() * "/../../../data/breast_21/originalGenomes.txt", ',', Float64, '\n')

#### nmf algorithm
nmf_global_params = NMFParams(
    init          = :random,
    dist          = Uniform(0, 1),
    max_iter      = 10_000,
    stopping_crit = :conn,
    stopping_iter = 100,
    verbose       = true,
)

# nmf_local_params = MUParams(
#     β       = 1,
#     div     = :β,
#     scale_W = false,
#     alg     = :mu
# )

nmf_local_params = FIParams(
    β       = 1,
    scale_W = false,
    alg     = :mm
)

W_fit_julia = zeros(Float64, size(W_fit))
H_fit_julia = zeros(Float64, size(H_fit))

for j =  1:10
    V_iter = readdlm(pwd() * "/../../../data/breast_21/V$(j)_rank$(K).txt", ',', Float64, '\n')

    w_init = reshape(W_init[:,j], size(W_init,1), 1)
    h_init = reshape(H_init[j,:], 1, size(H_init,2))

    res = VNMF._nmf_FI(V_iter, w_init, h_init, nmf_global_params, nmf_local_params)

    VNMF.scale_col_W!(res.W, res.H)

    W_fit_julia[:, j] .= reshape(res.W, length(res.W))
    H_fit_julia[j, :] .= reshape(res.H, length(res.H))
end

@test maximum(abs.(W_fit_julia - W_fit)) ≈ 0 atol = 1e-10
@test maximum(abs.(H_fit_julia - H_fit)) ≈ 0 atol = 1e-10

##### ANALYSE W_ALL AND H_ALL
W_all = W_fit_julia
H_all = H_fit_julia

##### REPRODUCE FOR RANK = 5

nmf_global_params = NMFParams(
    init          = :random,
    dist          = Uniform(0, 1),
    rank          = 5,
    max_iter      = 20_000,
    stopping_crit = :cost,
    stopping_tol  = 1e-6,
    verbose       = false,
)

nmf_local_params = MUParams(
    β       = 2,
    div     = :β,
    scale_W = true,
    alg     = :mu
)

nmf = NMF(
    solver        = nmf_MU,
    global_params = nmf_global_params,
    local_params  = nmf_local_params
)

#### build grid
list_V_pert = VNMF._get_list_perturbed_matrix(V, 50, :Multinomial, rng=rng)
list_nmf    = Vector{NMF}()
for (rank, seed) in zip(5:5, 5:5)
    #### using a different seed for each run of NMF
    nmf_rank = copy(nmf)
    nmf_rank.global_params.rank = rank
    nmf_rank.global_params.seed = seed

    push!(list_nmf, nmf_rank)
end

grid = vec(collect(Iterators.product(list_V_pert, list_nmf)))

#### run in parallel
grid_nmf_results = @showprogress pmap(grid) do g
    VNMF._rank_by_stability_one(g)
end

W_all = hcat([nmf_results.W for nmf_results in grid_nmf_results]...)
H_all = vcat([nmf_results.H for nmf_results in grid_nmf_results]...)

#### recover parameters of the algorithm
n_fac = nmf_global_params.rank
n_nmf = size(grid_nmf_results, 1)
n_tot = n_fac * n_nmf


text:
# Quadratic programming in MATLAB-like style
#
#   minimize x^2 + xy + y^2 + yz + z^2
#
#   s.t.    x + 2 y + 3 z >= 4
#           x +   y       >= 1
#

using Gurobi, Test

@testset "QP2" begin

	env = Gurobi.Env()
	setparam!(env, "OutputFlag", 0)

	model = gurobi_model(env;
		name = "qp_02",
		f = [0., 0., 0.],
		H = [2. 1. 0.; 1. 2. 1.; 0. 1. 2.],
		A = -[1. 2. 3.; 1. 1. 0.],
		b = -[4., 1.])

	optimize(model)

	@test isapprox(get_solution(model), [0.571429,0.428571,0.857143], atol=1e-5)
	@test isapprox(get_objval(model), 1.857143, atol=1e-5)
end


text:
"""
    linpred!(out, p::LinPred, f::Real=1.0)

Overwrite `out` with the linear predictor from `p` with factor `f`

The effective coefficient vector, `p.scratchbeta`, is evaluated as `p.beta0 .+ f * p.delbeta`,
and `out` is updated to `p.X * p.scratchbeta`
"""
function linpred!(out, p::LinPred, f::Real=1.)
    mul!(out, p.X, iszero(f) ? p.beta0 : broadcast!(muladd, p.scratchbeta, f, p.delbeta, p.beta0))
end

"""
    linpred(p::LinPred, f::Read=1.0)

Return the linear predictor `p.X * (p.beta0 .+ f * p.delbeta)`
"""
linpred(p::LinPred, f::Real=1.) = linpred!(Vector{eltype(p.X)}(undef, size(p.X, 1)), p, f)

"""
    installbeta!(p::LinPred, f::Real=1.0)

Install `pbeta0 .+= f * p.delbeta` and zero out `p.delbeta`.  Return the updated `p.beta0`.
"""
function installbeta!(p::LinPred, f::Real=1.)
    beta0 = p.beta0
    delbeta = p.delbeta
    @inbounds for i = eachindex(beta0,delbeta)
        beta0[i] += delbeta[i]*f
        delbeta[i] = 0
    end
    beta0
end

"""
    DensePredQR

A `LinPred` type with a dense, unpivoted QR decomposition of `X`

# Members

- `X`: Model matrix of size `n` × `p` with `n ≥ p`.  Should be full column rank.
- `beta0`: base coefficient vector of length `p`
- `delbeta`: increment to coefficient vector, also of length `p`
- `scratchbeta`: scratch vector of length `p`, used in `linpred!` method
- `qr`: a `QRCompactWY` object created from `X`, with optional row weights.
"""
mutable struct DensePredQR{T<:BlasReal} <: DensePred
    X::Matrix{T}                  # model matrix
    beta0::Vector{T}              # base coefficient vector
    delbeta::Vector{T}            # coefficient increment
    scratchbeta::Vector{T}
    qr::QRCompactWY{T}
    function DensePredQR{T}(X::Matrix{T}, beta0::Vector{T}) where T
        n, p = size(X)
        length(beta0) == p || throw(DimensionMismatch("length(β0) ≠ size(X,2)"))
        new{T}(X, beta0, zeros(T,p), zeros(T,p), qr(X))
    end
    function DensePredQR{T}(X::Matrix{T}) where T
        n, p = size(X)
        new{T}(X, zeros(T, p), zeros(T,p), zeros(T,p), qr(X))
    end
end
DensePredQR(X::Matrix, beta0::Vector) = DensePredQR{eltype(X)}(X, beta0)
DensePredQR(X::Matrix{T}) where T = DensePredQR{T}(X, zeros(T, size(X,2)))
convert(::Type{DensePredQR{T}}, X::Matrix{T}) where {T} = DensePredQR{T}(X, zeros(T, size(X, 2)))

"""
    delbeta!(p::LinPred, r::Vector)

Evaluate and return `p.delbeta` the increment to the coefficient vector from residual `r`
"""
function delbeta! end

function delbeta!(p::DensePredQR{T}, r::Vector{T}) where T<:BlasReal
    p.delbeta = p.qr\r
    return p
end

"""
    DensePredChol{T}

A `LinPred` type with a dense Cholesky factorization of `X'X`

# Members

- `X`: model matrix of size `n` × `p` with `n ≥ p`.  Should be full column rank.
- `beta0`: base coefficient vector of length `p`
- `delbeta`: increment to coefficient vector, also of length `p`
- `scratchbeta`: scratch vector of length `p`, used in `linpred!` method
- `chol`: a `Cholesky` object created from `X'X`, possibly using row weights.
- `scratchm1`: scratch Matrix{T} of the same size as `X`
- `scratchm2`: scratch Matrix{T} os the same size as `X'X`
"""
mutable struct DensePredChol{T<:BlasReal,C} <: DensePred
    X::Matrix{T}                   # model matrix
    beta0::Vector{T}               # base vector for coefficients
    delbeta::Vector{T}             # coefficient increment
    scratchbeta::Vector{T}
    chol::C
    scratchm1::Matrix{T}
    scratchm2::Matrix{T}
end
function DensePredChol(X::StridedMatrix, pivot::Bool)
    F = Hermitian(float(X'X))
    T = eltype(F)
    F = pivot ? cholesky!(F, Val(true), tol = -one(T), check = false) : cholesky!(F)
    DensePredChol(AbstractMatrix{T}(X),
        zeros(T, size(X, 2)),
        zeros(T, size(X, 2)),
        zeros(T, size(X, 2)),
        F,
        similar(X, T),
        similar(cholfactors(F)))
end

cholpred(X::StridedMatrix, pivot::Bool=false) = DensePredChol(X, pivot)

cholfactors(c::Union{Cholesky,CholeskyPivoted}) = c.factors
cholesky!(p::DensePredChol{T}) where {T<:FP} = p.chol

cholesky(p::DensePredQR{T}) where {T<:FP} = Cholesky{T,typeof(p.X)}(copy(p.qr.R), 'U', 0)
function cholesky(p::DensePredChol{T}) where T<:FP
    c = p.chol
    Cholesky(copy(cholfactors(c)), c.uplo, c.info)
end
cholesky!(p::DensePredQR{T}) where {T<:FP} = Cholesky{T,typeof(p.X)}(p.qr.R, 'U', 0)

function delbeta!(p::DensePredChol{T,<:Cholesky}, r::Vector{T}) where T<:BlasReal
    ldiv!(p.chol, mul!(p.delbeta, transpose(p.X), r))
    p
end

function delbeta!(p::DensePredChol{T,<:CholeskyPivoted}, r::Vector{T}) where T<:BlasReal
    ch = p.chol
    delbeta = mul!(p.delbeta, adjoint(p.X), r)
    rnk = rank(ch)
    if rnk == length(delbeta)
        ldiv!(ch, delbeta)
    else
        permute!(delbeta, ch.p)
        for k=(rnk+1):length(delbeta)
            delbeta[k] = -zero(T)
        end
        LAPACK.potrs!(ch.uplo, view(ch.factors, 1:rnk, 1:rnk), view(delbeta, 1:rnk))
        invpermute!(delbeta, ch.p)
    end
    p
end

function delbeta!(p::DensePredChol{T,<:Cholesky}, r::Vector{T}, wt::Vector{T}) where T<:BlasReal
    scr = mul!(p.scratchm1, Diagonal(wt), p.X)
    cholesky!(Hermitian(mul!(cholfactors(p.chol), transpose(scr), p.X), :U))
    mul!(p.delbeta, transpose(scr), r)
    ldiv!(p.chol, p.delbeta)
    p
end

function delbeta!(p::DensePredChol{T,<:CholeskyPivoted}, r::Vector{T}, wt::Vector{T}) where T<:BlasReal
    cf = cholfactors(p.chol)
    piv = p.chol.p
    cf .= mul!(p.scratchm2, adjoint(LinearAlgebra.mul!(p.scratchm1, Diagonal(wt), p.X)), p.X)[piv, piv]
    cholesky!(Hermitian(cf, Symbol(p.chol.uplo)))
    ldiv!(p.chol, mul!(p.delbeta, transpose(p.scratchm1), r))
    p
end

mutable struct SparsePredChol{T,M<:SparseMatrixCSC,C} <: GLM.LinPred
    X::M                           # model matrix
    Xt::M                          # X'
    beta0::Vector{T}               # base vector for coefficients
    delbeta::Vector{T}             # coefficient increment
    scratchbeta::Vector{T}
    chol::C
    scratch::M
end
function SparsePredChol(X::SparseMatrixCSC{T}) where T
    chol = cholesky(sparse(I, size(X, 2), size(X,2)))
    return SparsePredChol{eltype(X),typeof(X),typeof(chol)}(X,
        X',
        zeros(T, size(X, 2)),
        zeros(T, size(X, 2)),
        zeros(T, size(X, 2)),
        chol,
        similar(X))
end

cholpred(X::SparseMatrixCSC) = SparsePredChol(X)

function delbeta!(p::SparsePredChol{T}, r::Vector{T}, wt::Vector{T}) where T
    scr = mul!(p.scratch, Diagonal(wt), p.X)
    XtWX = p.Xt*scr
    c = p.chol = cholesky(Symmetric{eltype(XtWX),typeof(XtWX)}(XtWX, 'L'))
    p.delbeta = c \ mul!(p.delbeta, adjoint(scr), r)
end

LinearAlgebra.cholesky(p::SparsePredChol{T}) where {T} = copy(p.chol)
LinearAlgebra.cholesky!(p::SparsePredChol{T}) where {T} = p.chol

invchol(x::DensePred) = inv(cholesky!(x))
function invchol(x::DensePredChol{T,<: CholeskyPivoted}) where T
    ch = x.chol
    rnk = rank(ch)
    p = length(x.delbeta)
    rnk == p && return inv(ch)
    fac = ch.factors
    res = fill(convert(T, NaN), size(fac))
    for j in 1:rnk, i in 1:rnk
        res[i, j] = fac[i, j]
    end
    copytri!(LAPACK.potri!(ch.uplo, view(res, 1:rnk, 1:rnk)), ch.uplo, true)
    ipiv = invperm(ch.p)
    res[ipiv, ipiv]
end
invchol(x::SparsePredChol) = cholesky!(x) \ Matrix{Float64}(I, size(x.X, 2), size(x.X, 2))
vcov(x::LinPredModel) = rmul!(invchol(x.pp), dispersion(x, true))

function cor(x::LinPredModel)
    Σ = vcov(x)
    invstd = inv.(sqrt.(diag(Σ)))
    lmul!(Diagonal(invstd), rmul!(Σ, Diagonal(invstd)))
end

stderror(x::LinPredModel) = sqrt.(diag(vcov(x)))

function show(io::IO, obj::LinPredModel)
    println(io, "$(typeof(obj)):\n\nCoefficients:\n", coeftable(obj))
end

modelframe(obj::LinPredModel) = obj.fr
modelmatrix(obj::LinPredModel) = obj.pp.X
response(obj::LinPredModel) = obj.rr.y

fitted(m::LinPredModel) = m.rr.mu
predict(mm::LinPredModel) = fitted(mm)
formula(obj::LinPredModel) = modelframe(obj).formula
residuals(obj::LinPredModel) = residuals(obj.rr)

"""
    nobs(obj::LinearModel)
    nobs(obj::GLM)

For linear and generalized linear models, returns the number of rows, or,
when prior weights are specified, the sum of weights.
"""
function nobs(obj::LinPredModel)
    if isempty(obj.rr.wts)
        oftype(sum(one(eltype(obj.rr.wts))), length(obj.rr.y))
    else
        sum(obj.rr.wts)
    end
end

coef(x::LinPred) = x.beta0
coef(obj::LinPredModel) = coef(obj.pp)

dof_residual(obj::LinPredModel) = nobs(obj) - dof(obj) + 1


text:
using MCParticles
using Documenter

makedocs(;
    modules=[MCParticles],
    authors="Joshua Bon <joshuajbon@gmail.com>",
    repo="https://github.com/bonStats/MCParticles.jl/blob/{commit}{path}#L{line}",
    sitename="MCParticles.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://bonStats.github.io/MCParticles.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/bonStats/MCParticles.jl",
)


text:
# all things geodesic including transformations between standard
# geospatial coordinate systems
using Geodesy

# for direct manipulations of Cartesian and Longitude-Latitude coordinates
using CoordinateTransformations, Rotations, StaticArrays

# geometric and trigonomic constants
include("mathies.jl")

Φ
icoso_circumradius = sqrt(Φ^2 + 1)
r = earth_radius = 6378137.0

# pi/10 is 18 degrees, 3 * pi/10 is 54
# in geographic terms, the first element of each pair is North (Z+),
#   the second is East (Y+)
directions = SVector.([ ( 0           ,  0            ),
                        ( sin(pi/10)  , -cos(pi/10)   ),
                        (-sin(3*pi/10), -cos(3*pi/10) ),
                        (-sin(3*pi/10),  cos(3*pi/10) ),
                        ( sin(pi/10)  ,  cos(pi/10)   ),
                        ( Φ           ,  0            )])
proportioned_dirs = r * directions

rotY = LinearMap(RotZ(directions[6][1]))
rotZ = LinearMap(RotY(directions[6][2]))
rotZY = compose(rotZ, rotY)
rotYZ = compose(rotY, rotZ)

protY = LinearMap(RotZ(proportioned_dirs[6][1]))
protZ = LinearMap(RotY(proportioned_dirs[6][2]))
protZY = compose(rotZ, rotY)
protYZ = compose(rotY, rotZ)

composite = LinearMap(RotZ(directions[2][2])) ∘ LinearMap(RotY(directions[2][1]))
comp2 = inv(composite)

old_ecef = ECEF(0.0, earth_radius * -icoso_circumradius, 0.0)

tolla = LLAfromECEF(wgs84)

old_lla = tolla(old_ecef)

rotz_old = tolla(ECEF(rotZ(old_ecef)))
roty_old = tolla(ECEF(rotY(old_ecef)))

rotzy_old = tolla(ECEF(rotZY(old_ecef)))
rotyz_old = tolla(ECEF(rotYZ(old_ecef)))

protz_old = tolla(ECEF(protZ(old_ecef)))
proty_old = tolla(ECEF(protY(old_ecef)))

protzy_old = tolla(ECEF(protZY(old_ecef)))
protyz_old = tolla(ECEF(protYZ(old_ecef)))

max_stars_repo_name: sespaces/modjin

max_issues_repo_name: sespaces/modjin

max_forks_repo_name: sespaces/modjin


text:
module DataEditing

    using DataFrames
    using CSV
    using DataFramesMeta
    using Query
    # local
    using ONSCodes

    # local - run add_scripts_to_path.jl first

    using GlobalDecls
    using Utils

    export make_edited_datasets, add_ons_codes_to!,make_edited_hb_data, get_ofsted_data
    export create_actual_2020_allocation, parse_2020_allocation

    function parse_2020_allocation( filename :: String ) :: Dict
        n = 0
        key = "x"
        out=Dict{String,NamedTuple}()
        oldv = -1
        newv = -1
        for line in eachline( filename )
            n += 1
            r = (n % 3)
            if r==1
                  key = line
                  if key == "#"
                      break;
                  end
            elseif r == 2
                oldv = parse(Int64, line )
            else
                newv = parse(Int64, line )
                out[key]=(oldv=oldv,newv=newv)
            end;
            print("$key\n");
        end # loop
        out
    end # parse_allocation

    function create_actual_2020_allocation()
        start_year = 2020
        dd = DATADIR*"edited/$(start_year)"
        gfname = "$dd/GRANTS_$(start_year).csv"
        grantdata = CSV.File( gfname ) |> DataFrame
        grantdata[!,:my_estimated_amount] .= 0
        grantdata[!,:my_estimated_amount]= copy(grantdata[!,:amount])
        new_grants = parse_2020_allocation( "$dd/raw_grant_data_tmp.txt" );
        for (k,v) in new_grants
            grantdata[(grantdata.council.==k),:amount] .= v.newv
        end
        CSV.write( gfname, grantdata )
        grantdata
    end

    function yearstr( year :: Integer ) :: String
        ys = "?"
        if year == 2018
            ys = "2017-18"
        elseif year == 2019
            ys = "2018-19"
        end
        @assert ys != "?"
        ys
    end

    """
    FIXME we only actually have 1 set of this data
    """
    function get_payment_data( year :: Integer )::DataFrame
        CSV.File(
            DATADIR*"Staying_Put_grant_letter_to_LAs_2019-20_220319.csv",
            delim=',' ) |> DataFrame
    end

    """
    This is the LA payment spreadsheet from [here]()
    """
    PRE_PAYMENT_DATA = get_payment_data( GlobalDecls.SIMULATION_YEAR )

    function get_ofsted_data( year :: Integer ) :: DataFrame
        ys = yearstr( year )
        CSV.File(
            DATADIR*"Fostering_in_England_$(ys)_dataset_transposed_extended.csv",
            delim=',',
            missingstrings=["-999","",".."],
            types=make_type_block(3:1000,Float64)
             ) |> DataFrame
    end

    """
    This is my augmented version of the OFSTED fostering spreadsheet from [here](), with added info
    from the [AFC Fostering Study Spreadsheet]().
    """
    PRE_OFDATA =  get_ofsted_data( GlobalDecls.SIMULATION_YEAR )

    function fixup_cc_codes()
        n = size(LAMAPPINGS)[1]
        fixes = DataFrame( CSV.File(
           DATADIR*"las/Local_Authority_Districts_April_2015_Names_and_Codes_in_the_United_Kingdom.csv" ))
        for i in 1:n
            mapname = LAMAPPINGS[:name][i];
            mapcode = LAMAPPINGS[:new_gss_code][i]
            if mapname in fixes[:LAD15NM]
               fixrec = fixes[(fixes.LAD15NM.==mapname),:];
               fixcode = fixrec[:LAD15CD][1];
               # println( "fixcode $fixcode mapcode $mapcode ")
               if fixcode != mapcode
                   println( "change $mapcode to $fixcode for $mapname" );LAMAPPINGS[:new_gss_code][i] = fixcode
               end
            end
        end
    end # fixup_cc_codes

    """
    add ONS codes for the LA and enclosing region to some dataframe with lanames
    * `cname` - fieldname with LA names in it
    * `ccode` - name for the added ONS LA code
    * `rcode` - name for ONS region code field

    returns - list of names that have not been mapped
    """
    function add_ons_codes_to!( data :: DataFrame, cname :: Symbol, ccode :: Symbol, rcode :: Symbol )
        println( data[cname] )
        data[ccode]=map( c->code_from_name(c), data[cname])
        data[rcode]=map( c->region_code_from_name(c), data[cname])
    end

    function make_edited_datasets()
        path = DATADIR*"edited/"
        add_ons_codes_to!( PRE_PAYMENT_DATA, :council, :ccode, :rcode )
        CSV.write( path*"PAYMENT_DATA.csv", PRE_PAYMENT_DATA )
        add_ons_codes_to!( PRE_OFDATA, :council, :ccode, :rcode )
        CSV.write( path*"OFDATA.csv", PRE_OFDATA )
    end

    function make_edited_hb_data( year )
        # just for reference ..
        if year == 2019
            hbdata =  load( DATADIR*"2019-20_LHA_TABLES.csv" ) |> DataFrame
            add_ons_codes_to!( hbdata, :BRMA, :ccode, :rcode )
            CSV.write( DATADIR*"edited/"*"LHA_2019.csv", hbdata )
        elseif year == 2020
            hbdata =  load( DATADIR*"202021_LHA_TABLES.csv" ) |> DataFrame
            add_ons_codes_to!( hbdata, :BRMA, :ccode, :rcode )
            CSV.write( DATADIR*"edited/"*"LHA_2020.csv", hbdata )
        end
    end

end # DataEditing


text:
function Dmat(N, m, lam)

    # Initialise D with identity or half-derivative accordingly:
    if (round(m) == m) && (lam == 1/2)
        e = (2^m*gamma(m+1/2)/sqrt(pi)) .* ones(N)
        D = myspdiagm(spdiagm(m => e), N)
    elseif (round(m) !== m) && (lam == 1)
        e = 2^(m-1/2)*gamma(m+1) * ones(N,1)
        D = myspdiagm(spdiagm(Int(m-1/2) => e, Int(m+1/2) =>e), N)
    elseif round(m) == m
        D = 1
        for k = 0:floor(m)-1
            D = Dmat1(N, -k+0.5, lam+k) * D
        end
    else
        D = Dmat05(N, lam)
        for k = 0:floor(m)-1
            D = Dmat1(N, -k-0.5, lam+k+0.5) * D
        end
    end
    return D
end

function Dmat05(N, lam)
    # Half derivative: d^{1/2}/dx^{1/2} (1+x)^(lam-1/2)*C^{(lam)}_n(x)
    e = (gamma(lam+0.5)/gamma(lam)) .* ones(N)
    D = myspdiagm(spdiagm(0 => e, 1 => e), N)
    return D
end

function Dmat1(N, mu, lam)
    # Full derivative: d/dx (1+x)^mu*C^{(l)}_n(x)
    e = ones(N)
    if mu == 0
        D = 2*lam*myspdiagm(spdiagm(1 => e), N)
    else
        v = (mu-lam)./(collect(0:N-1).+lam)
        D = myspdiagm(spdiagm(0 => lam.*(1 .+v), 1 => lam.*2 .*e, 2 => (lam.*(1 .-v))[3:end]), N)
    end
end

myspdiagm(M, N) = size(M, 1) !== N ? (return M[1:N, 1:N]) : (return M)


text:
# module ODESolversTests

using GridapODEs
using GridapODEs.ODETools: GenericODESolution
using GridapODEs.ODETools: BackwardEuler
using GridapODEs.ODETools: RungeKutta
using GridapODEs.ODETools: ThetaMethodNonlinearOperator
using GridapODEs.ODETools: solve!
using GridapODEs
using GridapODEs.ODETools
using Gridap
using Test

# using Gridap.Algebra: residual, jacobian

include("ODEOperatorMocks.jl")

op = ODEOperatorMock{Float64,Constant}(1.0,0.0,1.0)

include("ODESolverMocks.jl")

t0 = 0.0
tf = 1.0
dt = 0.1
u0 = ones(2)*2

# NonlinearOperator tests

sop = OperatorMock(op,tf,dt,u0)
isa(sop,NonlinearOperator)

ode_cache = allocate_cache(op)

x = zero_initial_guess(sop)
x .+= 1.0
isa(sop,OperatorMock)
isa(x,AbstractVector)
r = allocate_residual(sop,x)
J = allocate_jacobian(sop,x)
residual!(r,sop,x)
jacobian!(J,sop,x)
@test all(r .== [ -11.0 -11.0])
@test all(J .== [ 9.0 0.0; 0.0 9.0])
_r = residual(sop,x)
_J = jacobian(sop,x)
@test all(_r .== [ -11.0 -11.0])
@test all(_J .== [ 9.0 0.0; 0.0 9.0])

# NLSolver tests

nls = NLSolverMock()
cache = solve!(x,nls,sop)
r, J, dx = cache
@test all(r.==_r)
@test all(J.==_J)
@test all(dx.≈11/9)
@test all(x.≈1+11/9)

#ODESolver tests

odesol = ODESolverMock(nls,dt)
uf = copy(u0)
uf.=1.0

uf, tf, cache = solve_step!(uf,odesol,op,u0,t0,nothing)
uf
@test tf==t0+dt
@test all(uf.≈x)

# ODESolutions

tF = 10.0
sol = GenericODESolution(odesol,op,u0,t0,tF)
current, state = Base.iterate(sol)
uf, tf = current
@test tf==t0+dt
@test all(uf.≈x)

# BackwardEulerNonlinearOperator tests

tf = t0+dt
vf = copy(u0)
sop = ThetaMethodNonlinearOperator(op,tf,dt,u0,ode_cache,vf) # See below
x = zero_initial_guess(sop)
x .+= 1.0
r = allocate_residual(sop,x)
J = allocate_jacobian(sop,x)
residual!(r,sop,x)
jacobian!(J,sop,x)
@test all(r .== [ -11.0 -11.0])
@test all(J .== [ 9.0 0.0; 0.0 9.0])
_r = residual(sop,x)
_J = jacobian(sop,x)
@test all(_r .== [ -11.0 -11.0])
@test all(_J .== [ 9.0 0.0; 0.0 9.0])

# BackwardEuler tests

ls = LUSolver()
odesol = BackwardEuler(ls,dt)
uf = copy(u0)
uf.=1.0
cache = nothing
# Juno.@enter solve_step!(uf,odesol,op,u0,t0,cache)
uf, tf, cache = solve_step!(uf,odesol,op,u0,t0,cache)
@test tf==t0+dt
@test all(uf.≈1+11/9)

@test test_ode_solver(odesol,op,u0,t0,tf)
test_ode_solver(odesol,op,u0,t0,tf)

# Affine and nonlinear solvers

op = ODEOperatorMock{Float64,Nonlinear}(1.0,0.0,1.0)
cache = nothing
uf, tf, cache = solve_step!(uf,odesol,op,u0,t0,cache)
@test tf==t0+dt
@test all(uf.≈1+11/9)

op = ODEOperatorMock{Float64,Affine}(1.0,0.0,1.0)
cache = nothing
uf, tf, cache = solve_step!(uf,odesol,op,u0,t0,cache)
@test tf==t0+dt
@test all(uf.≈1+11/9)

# RK tests

ls = LUSolver()
# BE equivalent
odesol = RungeKutta(ls,dt,:BE_1_0_1)
uf = copy(u0)
uf.=1.0
cache = nothing
uf, tf, cache = solve_step!(uf,odesol,op,u0,t0,cache)
@test tf==t0+dt
@test all(uf.≈1+11/9)
# SDIRK 2nd order
odesol = RungeKutta(ls,dt,:SDIRK_2_1_2)
uf = copy(u0)
uf.=1.0
cache = nothing
uf, tf, cache = solve_step!(uf,odesol,op,u0,t0,cache)
@test tf==t0+dt
@test all(uf.≈u0*(1.0+dt/(2*(1-dt))+dt*(1-2*dt)/(2*(1-dt)^2)))
# TRBDF (2nd order with some 0 on the diagonal)
odesol = RungeKutta(ls,dt,:TRBDF2_3_3_2)
uf.=1.0
cache = nothing
uf, tf, cache = solve_step!(uf,odesol,op,u0,t0,cache)
@test tf==t0+dt
@test all(uf.≈u0*1.105215241)

@test test_ode_solver(odesol,op,u0,t0,tf)
test_ode_solver(odesol,op,u0,t0,tf)

# Newmark test

op = ODEOperatorMock{Float64,Constant}(1.0,0.0,0.0)
ls = LUSolver()
γ = 0.5
β = 0.25
odesol = Newmark(ls,dt,γ,β)
uf = copy(u0)
v0 = ones(2)*(β*dt)
a0 = 0.0*ones(2)
uf.=1.0
cache = nothing
# Juno.@enter solve_step!(uf,odesol,op,u0,t0,cache)
uf, tf, vf, af, cache = solve_step!(uf,odesol,op,u0,v0,a0,t0,cache)
aᵧ = γ*af .+ (1-γ)*a0
aᵦ = 2*β*af .+ (1-2*β)*a0
@test tf==t0+dt
@test all(vf .≈ (v0 + dt*aᵧ))
@test all(uf .≈ (u0 + dt*v0 + 0.5*dt^2*aᵦ))

# end #module


text:
# Autogenerated wrapper script for libmseed_jll for x86_64-apple-darwin
export libmseed

JLLWrappers.@generate_wrapper_header("libmseed")
JLLWrappers.@declare_library_product(libmseed, "@rpath/libmseed.3.0.10.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libmseed,
        "lib/libmseed.3.0.10.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()