text:
# November, 2020
# Kevin Dorma
# module for common hydraulic calculations
# rev 0, December 2020
# rev 1, February 2021, added network hydraulics

# I need to make some decisions about how to interace with these functions
# I assume that I use dataframes to store information about sizing lines and general hydraulics
# then I will execute a function and it will return a result, but not change the raw information
# size lines
    # Input dataframe for describing the lines
    # Output line descriptor, pipe schedule, size criteria DP/100 and C, IDmm for each criteria and the chosen IDmm
# pick NPS will return the next size larger for NPS.This will be the chosen size.
# we will have utility functions for increasing or decreasing the NPS by 1.

module Hydraulics2

using DataFrames
using CSV
using ExcelFiles # more convenient than CSV for user input

export calcMoodyF
export addPipeProperties, addFluidProperties, getReynolds, getSegmentDP
export getLineSize, getListLargerNPS, getLargerNPS
# and the network hydraulic files
export indexLookup, doNetworkHydraulics, extractResults, combineResults

# these are the reference data.
# and this needs to be loaded correctly when we use the code in a module
include("npsList.jl")
include("schedList.jl")
include("pipeRoughnessList.jl")
include("fitting3K.jl")
include("pipeIDlist.jl")

#schedList = CSV.File("schedList.csv") |> DataFrame
#pipeRoughness = CSV.File("pipeRoughness.csv") |> DataFrame
#fitting3K = CSV.File("fitting3K.csv") |> DataFrame
#IDmm = CSV.File("pipeIDlist.csv") |> DataFrame
## not used
#pipeTable = CSV.read("pipeIDtable.csv");

function packageInfo()
    # return information about the package as a string
    return("Hydraulics package, Kevin Dorma. Written in Julia, rev 1 February 2021.")
end

# need function to append piping data
# later

function calcMoodyF(Reynolds,eD)
    # Moody friction factor, this is Halaand
    # eD is e/D relative roughness
    # Reynolds is Reynolds number
    invSqrtF = -1.8 .* log10.((eD ./ 3.70) .^ 1.11 + 6.9 ./ Reynolds)
    moodyF = (1 ./ (invSqrtF .^2) )
    return (moodyF)
end

function checkFittingList(fittingList)
    # fittingList is our list of fittings, fittingReference is the generic data
    # go through the list of fittings (fittingList) and compare with the reference list (possibly our fitting3K global list)
    # return items that are not found
    errorList = DataFrame(message = String[], entry=Int64[], Segment=String[], fittingType=String[])
    for i = 1:(size(fittingList)[1])
        thisFitting = fittingList[i,:fittingType]
        match = fitting3K[fitting3K.fittingType .== thisFitting,:]
        if ((size(match)[1]) == 0)
            push!(errorList, ("Fitting not found in row", i, fittingList[i,:Segment], fittingList[i,:fittingType] ))
        end
    end
    return (errorList)
end

function checkLineList(lines,fluidList)
    # go through our list of lines and compare with our reference lists
    # return the line items that are not found
    errorList = DataFrame(message = String[], entry=Int64[], Segment=String[], NPS=Float64[], Schedule=String[], material=String[], fluidName=String[])
    for i = 1:(size(lines)[1])
        # NPS, must check if this field exists
	if (indexLookup("NPS",names(lines)) > 0)
        thisItem = lines[i,:NPS]
        match = npsList[npsList.NPS .== thisItem,:]
        if ((size(match)[1]) == 0)
            push!(errorList, ("NPS not found in row", i, lines[i,:Segment], lines[i,:NPS], lines[i,:Schedule], lines[i,:Material], lines[i,:fluidName] ));
        end;
	end;
        # schedule
        thisItem = lines[i,:Schedule]
        match = schedList[schedList.Schedule .== thisItem,:]
        if ((size(match)[1]) == 0)
            push!(errorList, ("Schedule not found in row", i, lines[i,:Segment], lines[i,:NPS], lines[i,:Schedule], lines[i,:Material], lines[i,:fluidName] ))
        end
        # roughness material
        thisItem = lines[i,:Material]
        match = pipeRoughness[pipeRoughness.Material .== thisItem,:]
        if ((size(match)[1]) == 0)
            push!(errorList, ("Material not found in row", i, lines[i,:Segment], lines[i,:NPS], lines[i,:Schedule], lines[i,:Material], lines[i,:fluidName] ))
        end
        # fluid
        thisItem = lines[i,:fluidName]
        match = fluidList[fluidList.fluidName .== thisItem,:]
        if ((size(match)[1]) == 0)
            push!(errorList, ("FluidName not found in row", i, lines[i,:Segment], lines[i,:NPS], lines[i,:Schedule], lines[i,:Material], lines[i,:fluidName] ))
        end
    end
    return (errorList)
end

function addPipeProperties(df)
    # roughness, IDmm, eD ratio
    # this appends columns to df
    # used either line sizing or hydraulic calculations

    df[:,:roughnessMM] .= 0.0
    df[:,:IDmm] .= 0.0

    for i = 1:(size(df)[1])
        df[i,:roughnessMM]=pipeRoughness[df[i,:Material] .== pipeRoughness[:,:Material],:roughnessMM][1]
        df[i,:IDmm]= IDmm[(df[i,:NPS] .== IDmm[:,:NPS]) .& (df[i,:Schedule] .== IDmm[:,:Schedule]),:Idmm][1]
    end

    return (0.0)
end

function addFluidProperties(df,fluidList)
    # extract the density, viscosity and whatever else is needed from the fluidList
    # Add this to the df
    # used either line sizing or hydraulic calculations

    df[:,:rho_kgm3] .= 0.0
    df[:,:mu_mPas] .= 0.0
    df[:,:roughnessMM] .= 0.0

    for i = 1:(size(df)[1])
        df[i,:rho_kgm3] = fluidList[df[i,:fluidName] .== fluidList[:,:fluidName], :rho_kgm3][1]
        df[i,:mu_mPas]  = fluidList[df[i,:fluidName] .== fluidList[:,:fluidName], :mu_mPas][1]
        df[i,:roughnessMM]  = pipeRoughness[df[i,:Material] .== pipeRoughness[:,:Material], :roughnessMM][1]
    end

    return (0.0)
end

function lineSizeDP(df)
    # this is for sizing lines based on pressure drop kPa per 100 m, return the ID in mm
    # very convenient form straight out of Perrys handbook
    # uses friction factor correlation from Chen
    # Do not modify the original data in df

    g = 9.80665
    Sf = df.kPaPer100m*1000 ./ (df.rho_kgm3*g*100)
    q = df.massFlow .* df.margin ./(df.rho_kgm3*3600)
    eps = (df.roughnessMM/1000)
    kinVisc = (df.mu_mPas/1000) ./ df.rho_kgm3
    termA = ((eps.^5)*g .* Sf ./ (q.^2)).^0.25
    termB = ((kinVisc .^ 5)./((q.^3) .* Sf * g )).^0.2
    termC = 0.125*(termA .+ termB).^0.2
    returnValue = (1000*(termC .* (q.^2) ./ (g*Sf)).^0.2)
    return (returnValue)
end

function lineSizeErosion(df)
    # df is the dataframe with hydraulic information
    # this is the erosion C factor method
    # C = v/sqrt(rho), with v in m/s and rho in kg/m3
    # for most service, C = 120, this is similar to C = 100 for imperial units
    # return the ID in mm

    g = 9.80665
    maxVeloc = df.frictionCsi ./ sqrt.(df.rho_kgm3)
    q = df.massFlow .* df.margin ./(df.rho_kgm3*3600)
    area = q ./ maxVeloc
    returnValue = (sqrt.(4*area ./ pi )*1000.0)
    return (returnValue)
end

function getLineSize(df,fluidList)
    # from the calculated line size for the two different methods, determine the required line size

    addFluidProperties(df,fluidList)

    theSchedule = df[:,:Schedule]

    theSegment = df[:,:Segment]

    tempDP100 = lineSizeDP(df)
    tempErosion = lineSizeErosion(df)
    tempNeeded = tempErosion[:] + tempDP100[:]
    for j in 1:size(tempDP100)[1]
        tempNeeded[j] =max(tempDP100[j], tempErosion[j])
    end
    returnValue = DataFrame(Segment = theSegment, mmDP100 = tempDP100, mmErosion = tempErosion, mmNeeded = tempNeeded, Schedule = theSchedule)
    return (returnValue)
end

# I need a function to pick the next larger line size given the pipe schedule

function getLargerNPS(ourIDmm, ourSchedule)
    # from the common list of line sizes, pick the pipe size that matches the required schedule
    # and is the next larger that the required ID
    # this works for a single line size
    ourScheduleList = IDmm[IDmm[:,:Schedule] .== ourSchedule, :];
    refinedList = ourScheduleList[(ourScheduleList[:,:Idmm] .- ourIDmm) .> 0.0,:];
    theMinID = minimum(refinedList[:,:Idmm])
    ourNPS = refinedList[(refinedList[:,:Idmm] .== theMinID), :NPS]

    return (ourNPS)
end

function getListLargerNPS(lineSizeDF)
    # we will use the required line size and schedule and return the NPS and actual ID
    diamList = DataFrame(Segment=String[], NPS=Float64[], Schedule=String[], IDmm=Float64[])
    for i = 1:(size(lineSizeDF)[1])
        ourSchedule = lineSizeDF[i,:Schedule];
        ourIDmm = lineSizeDF[i,:mmNeeded]
        ourScheduleList = IDmm[IDmm[:,:Schedule] .== ourSchedule, :];
        refinedList = ourScheduleList[(ourScheduleList[:,:Idmm] .- ourIDmm) .> 0.0,:];
        theMinID = minimum(refinedList[:,:Idmm])
        ourNPS = refinedList[(refinedList[:,:Idmm] .== theMinID), :NPS][1]

        push!(diamList, (lineSizeDF[i,:Segment], ourNPS, ourSchedule, theMinID))
    end
    return (diamList)
end

function sizing2hydraulics(sizingDF, chosenLineSizeDF)
    # given the line sizing DF, copy the info into hydraulics dataframe
    # chosenSizeDF contains the NPS and Schedule
    # return the hydraulics DF
    # we also need a simple fittingDF, but this is done with a separate function
    hydraulicsDF = DataFrame(Segment=String[], Description=String[], LineTag=String[], PnID=String[], NPS=Float64[], Schedule=String[],
		Material=String[],	fluidName=[],	inletP_kPaa=[],	massFlow=Float64[])
    for i = 1:(size(sizingDF)[1])
        push!(hydraulicsDF, (sizingDF[i,:Segment], sizingDF[i,:Description],
                sizingDF[i,:LineTag], sizingDF[i,:PnID], chosenLineSizeDF[i,:NPS],
                chosenLineSizeDF[i,:Schedule], sizingDF[i,:Material], sizingDF[i,:fluidName],
                100.0, (sizingDF[i,:massFlow]*sizingDF[i,:margin])))
    end
    return(hydraulicsDF)
end

function sizing2fittingList(sizingDF, chosenLineSizeDF)
    # given the line sizing DF, copy the info into the fitting list
    # chosenSizeDF contains the NPS and Schedule
    # return the pre-populsted fittign list DF
    fittingListDF = DataFrame(Segment=String[], fittingType=String[], num_length_m=Float64[], comment=String[], revision=String[])

    for i = 1:(size(sizingDF)[1])
        push!(fittingListDF, (sizingDF[i,:Segment], "PIPE", 100.0, "from line sizing", "A"))
    end
    return(fittingListDF)
end

function incrementLineSize(rowNum, lineSizeDF, increment)
    # increment the specified row number needed line size by +1 or -1
    # not the cleanest function, but it works

    neededDF = DataFrame(Segment=String[], mmNeeded=Float64[], Schedule=String[])

    i = rowNum
    push!(neededDF, (lineSizeDF[i,:Segment], lineSizeDF[i,:IDmm] + increment, lineSizeDF[i,:Schedule]))

    incrementedLines = getListLargerNPS(neededDF)

    lineSizeDF[i,:NPS] = incrementedLines[1,:NPS]
    lineSizeDF[i,:IDmm] = incrementedLines[1,:IDmm]

    return (3.14)
end

function getReynolds(df)
    # preliminary calculations for pipe hydraulics
    # return the velocity, Reynolds, friction factor, density, diameterMM, diameterInch

    theSegment = df.Segment
    diamMM = df.IDmm
    diamInch = diamMM / 25.4
    rho = df.rho_kgm3
    volFlow = df.massFlow ./ df.rho_kgm3
    area = 0.25 * pi * (df.IDmm / 1000.0).^2
    velocity = volFlow ./ area / 3600
    Reynolds = df.rho_kgm3 .* velocity .* (df.IDmm / 1000.0) ./ (df.mu_mPas/1000.0)
    eD = df.roughnessMM ./ df.IDmm
    moodyF = calcMoodyF(Reynolds, eD)
    returnValue = DataFrame(Segment = theSegment, velocity_ms = velocity, rho_kgm3 = rho, Re = Reynolds, frictF = moodyF, IDmm = diamMM, IDinch = diamInch)
    return (returnValue)
end

function elementDP(lines, fittingList)
    # given the hydraulic elements in df (long list of pipe and fittings)
    # and the preliminary values in prelim
    # calculate the DP in each element
    # dp = rho.f. (L/D) v2/2
    # dp = rho K v2/2
    # where I need to calculate K = K1/Re + Kinf*(1 + Kd/Dinch^0.3)
    # can I get all of the pipe segments and use K = Kp * f L/D, where Kp = 1 for pipe
    # return the value for Kp (1 for pipe, 0 for fitting) and the pressure drop in the item
    returnValue = DataFrame(Segment = String[], Kp = Float64[], fittingK=Float64[], pipeK=Float64[], elementK = Float64[], DPkpa = Float64[])

    # start with getting Reynolds and other important things
    prelim = getReynolds(lines)

    for i = 1:(size(fittingList)[1])
        theSegment = fittingList[i,:Segment]

        # first we get all of the K values
        valK1 = fitting3K[fitting3K[:,:fittingType] .== fittingList[i,:fittingType],:K1][1]
        valKinf = fitting3K[fitting3K[:,:fittingType] .== fittingList[i,:fittingType],:Kinf][1]
        valKd = fitting3K[fitting3K[:,:fittingType] .== fittingList[i,:fittingType],:Kd][1]
        valKp = fitting3K[fitting3K[:,:fittingType] .== fittingList[i,:fittingType],:Kp][1]

        # now we get the hydraulic properties
        ff = prelim[prelim[:,:Segment] .== fittingList[i,:Segment],:frictF][1]
        rho = prelim[prelim[:,:Segment] .== fittingList[i,:Segment],:rho_kgm3][1]
        idInch = prelim[prelim[:,:Segment] .== fittingList[i,:Segment],:IDinch][1]
        idMM = prelim[prelim[:,:Segment] .== fittingList[i,:Segment],:IDmm][1]
        veloc = prelim[prelim[:,:Segment] .== fittingList[i,:Segment],:velocity_ms][1]
        Re = prelim[prelim[:,:Segment] .== fittingList[i,:Segment],:Re][1]

        Kfitting = valK1/Re + valKinf*(1.0 + valKd/(idInch^0.3))

        # pressure drop in kPa
        thisKfitting = fittingList[i,:num_length_m] * Kfitting
        thisKpipe = valKp * ff * fittingList[i,:num_length_m] /(idMM/1000.0)
        thisK = thisKfitting + thisKpipe
        thisDP = thisK * (0.001 * rho * 0.5 * (veloc)^ 2)

        push!(returnValue, (theSegment, valKp, thisKfitting, thisKpipe, thisK, thisDP))
    end
    return (returnValue)
end

function compileDP(lines, fittingDP)
    # for each entry in lines, sum all of the fittingDP
    returnValue = DataFrame(Segment = String[], segmentK = Float64[], DPkpa = Float64[], inletP = Float64[], outletP = Float64[])

    for i = 1:(size(lines)[1])
        theSegment = lines[i,:Segment]
        theSumK = sum(fittingDP[fittingDP[:,:Segment] .== theSegment,:elementK])
        theSumDP = sum(fittingDP[fittingDP[:,:Segment] .== theSegment,:DPkpa])
        theInletP = lines[i,:inletP_kPaa]
        theOutletP = theInletP - theSumDP
        push!(returnValue, (theSegment, theSumK, theSumDP, theInletP, theOutletP))
    end
    return (returnValue)
end

function getSegmentDP(lines, fittingList, fluidList)
    # lines describes the hydraulics, fittingList has the piping details, fluidList is the fluid

    addFluidProperties(lines,fluidList)
    addPipeProperties(lines)
    fittingDP = elementDP(lines, fittingList)
    returnValue = compileDP(lines, fittingDP)
    return (returnValue)
end

# functions added for network hydraulics
# february, 2021
# exported
function indexLookup(findMe, inThis)
    # find the index where findMe can be found inThis, which is a vector
    maxI = size(inThis)[1]
    foundIdx = 0
    for i = 1:maxI
        if (findMe == inThis[i])
            foundIdx = i
        end
    end
    return(foundIdx)
end

function getConnectivity(connectivity, nodeList)
    # connectivity, lineList and nodeList are dataframe
#    maxLine = size(lineList)[1]

    maxLine = size(connectivity)[1]

    maxNode = size(nodeList)[1]
    C = zeros(Float64,maxLine,maxNode)

    maxConnect = size(connectivity)[1]

    for i = 1:maxConnect
        theSegment = connectivity[i,:Segment]
        inNode = connectivity[i,:inNode]
        outNode = connectivity[i,:outNode]
#        theRow = indexLookup(theSegment,lineList[:,:Segment])
        theRow = indexLookup(theSegment,connectivity[:,:Segment])
        inCol = indexLookup(inNode,nodeList[:,:Node])
        outCol = indexLookup(outNode,nodeList[:,:Node])
        C[theRow,inCol] = 1.0
        C[theRow,outCol] = -1.0
    end
    return(C)
end

function removeTrivialRows(A)
    # remove rows from matrix A that contain only one entry
    numRows = size(A)[1]
    numColumns = size(A)[2]
    indexList = zeros(Int64, numRows,1)
    numNonTrivialRows = 0
    for i = 1:numRows
        countNonZero = 0
        for j = 1:numColumns
            if (A[i,j] != 0.0)
                countNonZero = countNonZero + 1
            end
        end
        if (countNonZero >= 2)
            numNonTrivialRows += 1
            indexList[numNonTrivialRows] = i
        end

    end
    returnMatrix = zeros(numNonTrivialRows,numColumns)
    for i = 1:numNonTrivialRows
        for j = 1:numColumns
            returnMatrix[i,j] = A[indexList[i],j]
        end
    end

    return(returnMatrix)
end

function initializeHydraulicParameters(setRe, df)
    # set all of the Reynolds numbers in lines df
    # return the hydraulicParameters (Reynolds, friction factor, density, diameterMM, diameterInch)

    theSegment = df.Segment
    diamMM = df.IDmm
    rho_kgm3 = df.rho_kgm3
    diamInch = diamMM / 25.4

    volFlow = df.massFlow ./ df.rho_kgm3
    area = 0.25 * pi * (df.IDmm / 1000.0).^2
    velocity = volFlow ./ area / 3600
    Reynolds = setRe .* (df.rho_kgm3 ./ df.rho_kgm3)
    eD = df.roughnessMM ./ df.IDmm
    moodyF = calcMoodyF(Reynolds, eD)

    returnValue = DataFrame(Segment = theSegment, rho_kgm3 = rho_kgm3, Re = Reynolds, frictF = moodyF, IDmm = diamMM, IDinch = diamInch)
    return (returnValue)
end

function initializeMassFlow(df, connectivity)
    # set mass flow in df to 1000 kg/h
    # but only if the flow rate is undefined or zero
    # I want to create a df with fields Segment and massFlow
    # we have viscosity, density and IDmm in the df
    # we will use a frictional pressure gradient of 20 kPa per 100 m for initial guess on mass flow
    # and friction factor 0.02
    # but then overwrite if we have something filled in for massFlow in the df

    returnValue = DataFrame(Segment = String[], massFlow = Float64[])

    # first, set all mass flows equal to 1000 kg/h
    numRows = size(connectivity)[1]
    for i = 1:numRows
        push!(returnValue, (connectivity[i,:Segment], 1000.0) )
    end

    # now, plow thorugh the hydraulic information
    # if there is no information here, then use DP per 100m to calculate a mass flow
    # if there is information here, use this
    numHydraulicRows = size(df)[1]

    dp100 = 20.0
    ff = 0.02 # estimate for friction factor

    for i = 1:numHydraulicRows
        ii = indexLookup(df[i,:Segment],connectivity[:,:Segment])
        # calculate this just in case
        dd = df[i,:IDmm]
        rho = df[i,:rho_kgm3]
        visc = df[i,:mu_mPas]
        mm = 3600 * sqrt(20*1000*rho*(pi^2 / (8*ff))*((dd^5)/100)*(1/1000)^5)
        if (ismissing(df[i,:massFlow]) || (df[i,:massFlow] == 0.0) )
            df[i,:massFlow] = mm
            returnValue[ii,:massFlow] = mm
        else
            returnValue[ii,:massFlow] = df[i,:massFlow]
        end
    end
    return (returnValue)
end

function getResistanceMatrix(hydraulicParameters, elementK, connectivity)
    # given the hydraulic parameters, a list of all segments
    # and the long list of K values for every element in the network
    # calculate the sum of the overall K for each segment
    # return a diagonal matrix with the K values on the diagonal
    # wait, the resistance K value is in terms of velocity in m/s
    # I need mass flow in kg/h
    #
    # I think I need to add the connectivity matrix to get the correct value for the variable number

    numRows = size(hydraulicParameters)[1]
    numColumns = size(connectivity)[1]

    numElements = size(elementK)[1]
#    matrixK = zeros(Float64,numRows,numRows)
#     this needs to be the number of variables in the connectivity matrix
    matrixK = zeros(Float64,numRows,numColumns)

    for i = 1:numElements
#        theRow = indexLookup(elementK[i,:Segment], hydraulicParameters[:,:Segment])
#        theRow = indexLookup(elementK[i,:Segment], hydraulicParameters[:,:Segment])

        theRow = indexLookup(elementK[i,:Segment], hydraulicParameters[:,:Segment])
        theColumn = indexLookup(elementK[i,:Segment], connectivity[:,:Segment])

        matrixK[theRow,theColumn] += elementK[i,:elementK]
    end
    # now we multiply each row element by 8 x 1000 / (rho pi^2 D^4) (1000/3600)^2
    # where D is in mm, mass flow in kg/h and DP in kPa
    for i = 1:numRows
        multK = 8.0*1000.0 * ((1000/3600)^2) / ((pi^2) * hydraulicParameters[i,:rho_kgm3] * (hydraulicParameters[i,:IDmm]^4)   )
        for j = 1:numColumns
            matrixK[i,j] *= multK
        end
    end
    return (matrixK)
end

function getSegLength(connectivity, fittingList)

    returnValue = DataFrame(Segment = String[], pipeLength=Float64[])

# initialize the segLength data
    for i=1:(size(connectivity)[1])
	push!(returnValue, (connectivity[i,:Segment], 0.0))
    end
    for i = 1:(size(fittingList)[1])
        valKp = fitting3K[fitting3K[:,:fittingType] .== fittingList[i,:fittingType],:Kp][1]
	ii = indexLookup(fittingList[i,:Segment], connectivity[:,:Segment])
	if (valKp == 1)
		returnValue[ii,:pipeLength] += fittingList[i,:num_length_m]
	end
    end
    return(returnValue)
end

function getConstraintMatrix(constraintDF, constraintType, segmentDF, nodeDF)
    # given the constraint data, the type (massFlow or pressure or value), the list of segments and nodes
    # create a matrix with values of 1 in the row, or the actual value

    numConstraints = size(constraintDF)[1]
    numSegments = size(segmentDF)[1]
    numNodes = size(nodeDF)[1]

    # we need the dimensions for the return matrix
    # assume a value for now
    numCols = 2
    if (constraintType == "flowRate")
        numCols = numSegments
    end
    if (constraintType == "pressure")
        numCols = numNodes
    end
    if (constraintType == "value")
        numCols = 1
    end
    returnMatrix = zeros(numConstraints, numCols)

    # cycle through all of the constraints
    # of the constraintType matches the :variable in the constraintDF, then place a 1 in the matrix entry
    for i = 1:numConstraints
        if (constraintType == "value")
            returnMatrix[i] = constraintDF[i,:value]
        elseif (constraintType == "flowRate")
            varIndexLocal = indexLookup(constraintDF[i,:object], segmentDF[:,:Segment])
            if (varIndexLocal > 0)
                returnMatrix[i,varIndexLocal] = 1.0
            end
        elseif (constraintType == "pressure")
            varIndexLocal = indexLookup(constraintDF[i,:object], nodeDF[:,:Node])
            if (varIndexLocal > 0)
                returnMatrix[i,varIndexLocal] = 1.0
            end
        end
    end

    return(returnMatrix)
end

function getQuadraticMatrix(hydraulicDF, which)
    # given the estimates of mass flow rate in the hydraulic matrix, fill out the
    # linearized estimate to the quadratic mass flow^2
    # m^2 ~= 2 M m - M*M,
    # or
    # m2 - 2 M m = - M*M, where we are looking for m2 (m^2) and m, and M is the
    # previous estiamte of m
    # the value of which is either quadratic, linear, value

    # given the constraint data, the type (massFlow or pressure or value), the list of segments and nodes
    # create a matrix with values of 1 in the row, or the actual value

    numEquations = size(hydraulicDF)[1]
    numCols = numEquations

    if (which == "rhs")
        numCols = 1
    end

    returnMatrix = zeros(numEquations, numCols)

    # cycle through all of the segments
    # fill in the value for the matrix equation
    # m2 - 2 M m = - M*M, where we are looking for m2 (m^2) and m, and M is the
    # if quadratic rM[i,i] = 1.0
    # if linear rm[i,i] = -2*hydraulicDF[i,:massFlow]
    # if value rm[i] = -hydraulicDF[i,:massFlow]*hydraulicDF[i,:massFlow]

    for i = 1:numEquations
        if (which == "rhs")
            returnMatrix[i] = -hydraulicDF[i,:massFlow]*hydraulicDF[i,:massFlow]
        elseif (which == "quadratic")
            returnMatrix[i,i] = 1.0
        elseif (which == "linear")
            returnMatrix[i,i] = -2.0 * hydraulicDF[i,:massFlow]
        end
    end

    return(returnMatrix)
end

function getSegmentIMatrix(hydraulicParameters, connectivity)
    # given the Identity matrix for the segments
    # return Identity matrix

    numRows = size(hydraulicParameters)[1]
    numColumns = size(connectivity)[1]

    matrixI = zeros(Float64,numRows,numColumns)

    for i = 1:numRows
        findMe = hydraulicParameters[i,:Segment]
        inThis = connectivity[:,:Segment]
        j = indexLookup(findMe, inThis)
        matrixI[i,j] = 1.0
    end

    return (matrixI)
end

function getElementK(prelim, fittingList)
    # given the hydraulic elements in df (long list of pipe and fittings)
    # and the preliminary friction properties (Re, eD, ff)
    # calculate the loss coefficient K in each element
    # where I need to calculate K = K1/Re + Kinf*(1 + Kd/Dinch^0.3)
    # can I get all of the pipe segments and use K = Kp * f L/D, where Kp = 1 for pipe
    # return the value for Kp (1 for pipe, 0 for fitting) and the pressure drop in the item
    returnValue = DataFrame(Segment = String[], Kp = Float64[], fittingK=Float64[], pipeK=Float64[], elementK = Float64[], pipeLength=Float64[])

    for i = 1:(size(fittingList)[1])
        theSegment = fittingList[i,:Segment]

        # first we get all of the K values
        valK1 = Hydraulics2.fitting3K[Hydraulics2.fitting3K[:,:fittingType] .== fittingList[i,:fittingType],:K1][1]
        valKinf = Hydraulics2.fitting3K[Hydraulics2.fitting3K[:,:fittingType] .== fittingList[i,:fittingType],:Kinf][1]
        valKd = Hydraulics2.fitting3K[Hydraulics2.fitting3K[:,:fittingType] .== fittingList[i,:fittingType],:Kd][1]
        valKp = Hydraulics2.fitting3K[Hydraulics2.fitting3K[:,:fittingType] .== fittingList[i,:fittingType],:Kp][1]

        # now we get the hydraulic properties
        ff = prelim[prelim[:,:Segment] .== fittingList[i,:Segment],:frictF][1]
        rho = prelim[prelim[:,:Segment] .== fittingList[i,:Segment],:rho_kgm3][1]
        idInch = prelim[prelim[:,:Segment] .== fittingList[i,:Segment],:IDinch][1]
        idMM = prelim[prelim[:,:Segment] .== fittingList[i,:Segment],:IDmm][1]
#        veloc = prelim[prelim[:,:Segment] .== fittingList[i,:Segment],:velocity_ms][1]
        Re = prelim[prelim[:,:Segment] .== fittingList[i,:Segment],:Re][1]

        Kfitting = valK1/Re + valKinf*(1.0 + valKd/(idInch^0.3))

        # pressure drop in kPa
        thisKfitting = fittingList[i,:num_length_m] * Kfitting
        thisKpipe = valKp * ff * fittingList[i,:num_length_m] /(idMM/1000.0)
	if (valKp == 1)
		thisPipeLength = fittingList[i,:num_length_m]
	else
		thisPipeLength = 0.0
	end

        thisK = thisKfitting + thisKpipe
#        thisDP = thisK * (0.001 * rho * 0.5 * (veloc)^ 2)

        push!(returnValue, (theSegment, valKp, thisKfitting, thisKpipe, thisK, thisPipeLength))
    end
    return (returnValue)
end

function addToBlockMatrix(gBM, localMatrix, gR, gC, eqnType, varType, numLocalRows, numLocalColumns)
    # we will traverse the localMatrix and place the entries in the global block matrix
    # gR and gC are the table of indices that map gR[eqnType,i] to the row in global matrix
    # and gC[varType,j] to the column in the global matrix
    # eqnType and varType refer to which block
    # numLocalRows and Columns tells us how many local entries are avalable

    for i = 1:numLocalRows
        for j = 1:numLocalColumns
            gi = gR[eqnType,i]
            gj = gC[varType,j]
            gBM[gi,gj] = localMatrix[i,j]
        end
    end
    return(3.14)
end

function addToBlockRHS(gBv, localVec, gR, eqnType, numLocalRows)
    # we will traverse the localMatrix and place the entries in the global block matrix
    # gR and gC are the table of indices that map gR[eqnType,i] to the row in global matrix
    # and gC[varType,j] to the column in the global matrix
    # eqnType and varType refer to which block
    # numLocalRows and Columns tells us how many local entries are avalable

    for i = 1:numLocalRows
        gi = gR[eqnType,i]
        gBv[gi] = localVec[i]
    end
    return(3.14)
end

# create routines to extract results
# extractFromBlockSolution(testSolution, globalColumn, j (which type is it),columnsInSegment[j])
function extractFromBlockSoln(gBx, gR, varType, numLocalRows)
    # we will traverse the appropriate entries in the global block matrix
    # gR is the vector of indices that map gR[varType,i] to the row in global solution vector
    # varType refer to which block
    # numLocalRows tells us how many local entries are avalable

    localX=zeros(numLocalRows,1)
    for i = 1:numLocalRows
        gi = gR[varType,i]
        localX[i] = gBx[gi]
    end
    return(localX)
end

# create routines to extract results
# extractFromBlockSolution(testSolution, globalColumn, j (which type is it),columnsInSegment[j])
function catIt(dfHydraulic1, dfHydraulic2, dfAll)
    # dfAll has all of the rows, we need to fill in blank rows for the other dfs

    localX=zeros(numLocalRows,1)
    for i = 1:numLocalRows
        gi = gR[varType,i]
        localX[i] = gBx[gi]
    end
    return(localX)
end

# create routines to extract results
function extractResults(gBx, lines, connectivity, nodeList, fittingList)
    # we will traverse the solutions in the block x vector
    # lines is the df with the segments defined
    # connectivity is the connectivity df
    # we will find all of the from entries in the connectivity matrix
    numLines = getSize(connectivity)
    numHydraulicLines = getSize(connectivity)
    numNodes = getSize(nodeList)

    columnsInSegment = getColumnsInSegment(numLines,numNodes)
    globalColumn = variableMapping(numLines, numNodes, columnsInSegment)

    ourDetails = (Hydraulics2.getReynolds(lines))
	# note that fittingList is a really long vector with every fitting
    ourElementK = getElementK(ourDetails, fittingList)
    Km = getResistanceMatrix(ourDetails, ourElementK, connectivity)  # this is an actual matrix
    segLength = getSegLength(connectivity,fittingList)  # dataframe for length and other stuff

    segList = copy(connectivity[:,:Segment])
    nList = copy(nodeList[:,:Node])

    # connectivity, used for pressure continuity
    localA = getConnectivity(connectivity,nodeList)

    # extract all of the entries with +1 and -1
    fromA = max.(localA, 0)
    toA = -min.(localA, 0)

    varType = 1        # mass flow
    xMassFlow = extractFromBlockSoln(gBx, globalColumn, varType, columnsInSegment[varType])

    varType = 2     # mass flow^2
    xMassFlow2 = extractFromBlockSoln(gBx,globalColumn,varType,columnsInSegment[varType])

    varType = 3     # DP
    xDP = extractFromBlockSoln(gBx,globalColumn,varType,columnsInSegment[varType])

    varType = 4     # node pressure
    xPnode = extractFromBlockSoln(gBx,globalColumn,varType,columnsInSegment[varType])
    fromP = fromA*xPnode
    toP = toA*xPnode

    # messy, but it works
#    allResults = [segList xMassFlow xDP fromP toP Konly]
    allResults = [segList xMassFlow xDP fromP toP]

#    returnVal = DataFrame(Segment=allResults[:,1], massFlow=allResults[:,2], dp=allResults[:,3], fromP=allResults[:,4], toP=allResults[:,5], elementK=allResults[:,6])

    returnVal = DataFrame(Segment=allResults[:,1], massFlow=allResults[:,2], dp=allResults[:,3], fromP=allResults[:,4], toP=allResults[:,5], fromNode=connectivity[:,:inNode], toNode=connectivity[:,:outNode], segmentLength=segLength[:,:pipeLength])

    return(returnVal)
end

function getSize(df)
    # get the number of rows
    return(size(df)[1])
end

function getColumnsInSegment(numLines, numNodes)
    # define the number of variables in each of the different types of variables
    numEqnTypes = 5 # not used
    numVarTypes = 4 # m, m2, dp, p

    columnsInSegment = zeros(Int64,numVarTypes,1)
    columnsInSegment[1] = numLines
    columnsInSegment[2] = numLines
    columnsInSegment[3] = numLines
    columnsInSegment[4] = numNodes

    return(columnsInSegment)
end

function variableMapping(numLines, numNodes, columnsInSegment)
    # create the variable mapping as a 2 D array
    # map[row,column] = global index
    # row is the type of variable (m, m2, dp, p)
    # column contains the index for each of the available variables in the global vector
    # function to count number of equations and variables
    numEqnTypes = 5 # not used
    numVarTypes = 4 # m, m2, dp, p

    maxOtherIndex = numLines + numNodes
    globalColumn = zeros(Int64,numVarTypes,maxOtherIndex)

    index = 1
    for i = 1:size(columnsInSegment)[1]
        for j = 1:columnsInSegment[i]
            globalColumn[i,j] = index
            index += 1
        end
    end

    return(globalColumn)
end

function doNetworkHydraulics(lineHydraulics, fluidList, connectivity, dofList, nodeList, fittingList)
    # there is a lot of code here. I do not want the user typing this. It should only be function calls
    #
    # check input data
        # numLines = checkSegments(lineHydraulics)
        # numNodes = checkNodes
        # I have a DOF check in the code
        # or numLines = getSize(lineHydraulics)
        # or numNodes = getSize(nodeList)
        # variableMapping = getVariableMapping(numLines, numNodes)
    #
    # resultVector = doHydraulics(lineHydraulics, fluidList,dofList,nodeList)
    # Then we need to create a new table with the segment results
    # and we need another table with the pressures
    #
    # then do this
    # segmentResults = extractResults("segment", lineHydraulics, resultVector)
    # nodeResults = extractResults("node", lineHydraulics, resultVector)
    # I could supply a list of dataframe names to include in the results
    #

    # driver script
    addFluidProperties(lineHydraulics,fluidList)
    addPipeProperties(lineHydraulics)

    # and this is where we need to initialize
    # return the matrix of mass flow and m2
    # I will need to provide the connectivity df as well to get all of the fields
    massFlow = initializeMassFlow(lineHydraulics,connectivity) # if we do not have a mass flow, make one up

    # pass the mass flow matrix
    #and set the Reynolds number and friction factor to a reasonable value
    hydraulicParameters = initializeHydraulicParameters(40000.0, lineHydraulics)

    numLines = size(connectivity)[1]
    numResistance = size(lineHydraulics)[1]
    numNodes = size(nodeList)[1];
    #C = zeros(Float64,numLines,numNodes)

    segList = copy(connectivity[:,:Segment])
#    segList = copy(lineHydraulics[:,:Segment])
    nList = copy(nodeList[:,:Node])

    # connectivity, used for pressure continuity
    # each row defines a segment: we have the inlet and outlet node
    # I don't think I need the lineHydraulics df
    Amatrix = getConnectivity(connectivity,nodeList)
    # I think I only need to assemble the block matrix.
    Ia = zeros(Float64,numLines,numLines)
    for i = 1:numLines
        Ia[i,i] = 1.0;
    end;

    # create the mass balance matrix
    Bprelim = Amatrix'
    Bmatrix = removeTrivialRows(Bprelim);

    # I need to create the constraint equations
    # use the connectivity df instead of lineHydraulics
    Cm = getConstraintMatrix(dofList, "flowRate", connectivity, nodeList)
    Cp = getConstraintMatrix(dofList, "pressure", connectivity, nodeList)
    xc = getConstraintMatrix(dofList, "value", connectivity, nodeList);

    # need to define where each matrix will go.
    # sure I could move this outside the loop, but it should not slow things down much
    # need vectors for the number of rows and columns in each segment
    # I could put most of this batch of code in a function
    # let
    # l = num lines
    # n = num nodes
    # k = resistance calcs
    # d = degrees of freedom

    # num unknonws = 3*l + n
    # count equations
    # num of mass balance relations = num rows in B
    # quadratic equations = l
    # resistance eqns = k
    # pressure connectivity = l
    # degree of freedom = d

        # the block matrix is
        #  1.  2.  3.  4.
        #    m, m2, dp,  p | rhs
        # 1  B,  0,  0,  0 | 0
        # 2 Qm,Qm2,  0,  0 | xq
        # 3  0,-Km,Iseg, 0 | 0
        # 4  0,  0, -I, Am | 0
        # 5 Cm,  0,  0, Cp | xc

    # function to count number of equations and variables
    # edit this to account fo rnumber of resistances
    numEqnTypes = 5
    numVarTypes = 4
    rowsInSegment = zeros(Int64,numEqnTypes,1)
    rowsInSegment[1] = size(Bmatrix)[1]
    rowsInSegment[2] = numLines
    rowsInSegment[3] = numResistance
    rowsInSegment[4] = numLines
    rowsInSegment[5] = size(dofList)[1]

    columnsInSegment = zeros(Int64,numVarTypes,1)
    columnsInSegment[1] = numLines
    columnsInSegment[2] = numLines
    columnsInSegment[3] = numLines
    columnsInSegment[4] = numNodes

    numEquations = sum(rowsInSegment)
    numVariables = sum(columnsInSegment)

    # I need this screening as a preliminary check
    (numEquations != numVariables) ? error("error: ", numEquations, " equations and ", numVariables, " variables") : println("OK: ", numEquations, " equations and variables")

    # this dozen lines could go into a function
    # I could create a matrix with global indices
    # globalRow[segment][localRow]
    # globalColumn[segment][localColumn]
    maxOtherIndex = numLines + numNodes
    globalRow = zeros(Int64,numEqnTypes,maxOtherIndex)
    globalColumn = zeros(Int64,numVarTypes,maxOtherIndex)

    index = 1
    for i = 1:size(rowsInSegment)[1]
        for j = 1:rowsInSegment[i]
            globalRow[i,j] = index
            index += 1
        end
    end

    index = 1
    for i = 1:size(columnsInSegment)[1]
        for j = 1:columnsInSegment[i]
            globalColumn[i,j] = index
            index += 1
        end
    end

    # define the global block matrix and vector

    nGlobalRow = maximum(globalRow)
    nGlobalColumn = maximum(globalColumn)
    gBM = zeros(Float64,nGlobalRow,nGlobalColumn) # global block matrix
    gBv = zeros(Float64,nGlobalRow,1) # global block vector
    gBx = zeros(Float64,nGlobalRow,1) # global block solution

    epsM = 0.01 # mass flow kg/h
    maxIter = 15

    errM = 100.0
    iter = 0

    # and this is where we need to iterate

    while ((errM > epsM) && (iter < maxIter) )
        # the resistance matrix is not applied to segments that are control valves
        ourElementK = getElementK(hydraulicParameters, fittingList)
        Km = getResistanceMatrix(hydraulicParameters, ourElementK, connectivity)
        Iseg = getSegmentIMatrix(hydraulicParameters, connectivity)

        # we need the quadratic matrix.
        # must realize that our lineHydraulic table is not the definative source of infomration
        # for this
        # we should use the connectivity table
        # Or do we create a simple table for only mass flow and m2?
        # m2 = 2 m M - M M
        # m2 - 2 m M = -M M
        # change this from lineHydraulics to massFlow
        # why does massFlow not work?
#        Qm2 = getQuadraticMatrix(lineHydraulics, "quadratic")
#        Qm  = getQuadraticMatrix(lineHydraulics, "linear")
#        xq  = getQuadraticMatrix(lineHydraulics, "rhs");

        Qm2 = getQuadraticMatrix(massFlow, "quadratic")
        Qm  = getQuadraticMatrix(massFlow, "linear")
        xq  = getQuadraticMatrix(massFlow, "rhs");

        # the block matrix is
        #  1.  2.  3.  4.
        #    m, m2, dp,  p | rhs
        # 1  B,  0,  0,  0 | 0
        # 2 Qm,Qm2,  0,  0 | xq
        # 3  0,-Km,Iseg, 0 | 0
        # 4  0,  0, -I, Am | 0
        # 5 Cm,  0,  0, Cp | xc
        #
        # P1 - P2 = K.m2, where P1 is in, P2 is out
        # DP = K.m2
        # P1 - P2 = DP
        # or
        # -K.m2 + DP = 0
        # -DP + P1 - P2 = 0

        # once I have the answer, how do I get inlet and outlet pressure for each segment.
            # The A matrix contains the information
            # The rows in A are for each segment.
            # columns with +1 are inlet, -1 are outlet.
            # Transform the A matrix by extracting the entries with +1. Then matrix multiply.
            # then get the entries with -1 and matrix multiply.

        # assemble the global block matrix
        i = 1
        j = 1
        addToBlockMatrix(gBM,Bmatrix,globalRow,globalColumn,i,j,rowsInSegment[i],columnsInSegment[j])

        i = 2
        j = 1
        addToBlockMatrix(gBM,Qm,globalRow,globalColumn,i,j,rowsInSegment[i],columnsInSegment[j])

        i = 2
        j = 2
        addToBlockMatrix(gBM,Qm2,globalRow,globalColumn,i,j,rowsInSegment[i],columnsInSegment[j])
        addToBlockRHS(gBv, xq, globalRow, i, rowsInSegment[i])

        i = 3
        j = 2
        addToBlockMatrix(gBM,-Km,globalRow,globalColumn,i,j,rowsInSegment[i],columnsInSegment[j])

        i = 3
        j = 3
        addToBlockMatrix(gBM,Iseg,globalRow,globalColumn,i,j,rowsInSegment[i],columnsInSegment[j])

        i = 4
        j = 3
        addToBlockMatrix(gBM,-Ia,globalRow,globalColumn,i,j,rowsInSegment[i],columnsInSegment[j])

        i = 4
        j = 4
        addToBlockMatrix(gBM,Amatrix,globalRow,globalColumn,i,j,rowsInSegment[i],columnsInSegment[j])

        i = 5
        j = 1
        addToBlockMatrix(gBM,Cm,globalRow,globalColumn,i,j,rowsInSegment[i],columnsInSegment[j])

        i = 5
        j = 4
        addToBlockMatrix(gBM,Cp,globalRow,globalColumn,i,j,rowsInSegment[i],columnsInSegment[j])
        addToBlockRHS(gBv, xc, globalRow, i, rowsInSegment[i])

        # and solve
        testSoln = gBM \ gBv

        # convergnece
#        errM = maximum(abs.(lineHydraulics[:,:massFlow] - testSoln[1:numLines]))

        errM = maximum(abs.(massFlow[:,:massFlow] - testSoln[1:numLines]))

        iter += 1

        # update the mass flow rate solutions
        massFlow[:,:massFlow] = testSoln[1:numLines]
        for i = 1:numResistance
            ii = indexLookup(lineHydraulics[i,:Segment], connectivity[:,:Segment])
            lineHydraulics[i,:massFlow] = testSoln[ii]
        end

        # then I need to update the hydraulic parameters
        hydraulicParameters = Hydraulics2.getReynolds(lineHydraulics)
        # and I need to update the pressures and DP estimates

        gBx = testSoln
    end
    return (gBx)
end

# create routines to extract results
# extractFromBlockSolution(testSolution, globalColumn, j (which type is it),columnsInSegment[j])
function combineResults(dfAllResults, dfLines)
    # dfAll has all of the rows
    # dfLines has the pressur drop hydraulic info
    # do an innerjoin to produce a single table

    hydraulicDetails = (Hydraulics2.getReynolds(dfLines));
    combined1=outerjoin(dfAllResults, hydraulicDetails, on = intersect(names(dfAllResults), names(hydraulicDetails)), makeunique=false,indicator=nothing,validate=(false,false));
    combined2=outerjoin(combined1, dfLines, on = intersect(names(combined1), names(dfLines)), makeunique=false,indicator=nothing,validate=(false,false));

    return(combined2)
end

end # module


text:
using Diagonalizations, LinearAlgebra, PosDefManifold, Test

# Method (1) real
n, t=10, 100
X=genDataMatrix(n, t)
Y=genDataMatrix(n, t)
Cx=Symmetric((X*X')/t)
Cy=Symmetric((Y*Y')/t)
Cxy=(X*Y')/t
cC=cca(Cx, Cy, Cxy, simple=true)
@test cC.F[1]'*Cx*cC.F[1]≈I
@test cC.F[2]'*Cy*cC.F[2]≈I
D=cC.F[1]'*Cxy*cC.F[2]
@test norm(D-Diagonal(D))+1. ≈ 1.

# Method (1) complex
Xc=genDataMatrix(ComplexF64, n, t)
Yc=genDataMatrix(ComplexF64, n, t)
Cxc=Hermitian((Xc*Xc')/t)
Cyc=Hermitian((Yc*Yc')/t)
Cxyc=(Xc*Yc')/t
cCc=cca(Cxc, Cyc, Cxyc, simple=true)
@test cCc.F[1]'*Cxc*cCc.F[1]≈I
@test cCc.F[2]'*Cyc*cCc.F[2]≈I
Dc=cCc.F[1]'*Cxyc*cCc.F[2]
@test norm(Dc-Diagonal(Dc))+1. ≈ 1.

# Method (2) real
cXY=cca(X, Y, simple=true)
@test cXY.F[1]'*Cx*cXY.F[1]≈I
@test cXY.F[2]'*Cy*cXY.F[2]≈I
D=cXY.F[1]'*Cxy*cXY.F[2]
@test norm(D-Diagonal(D))+1. ≈ 1.
@test cXY==cC

# Method (2) complex
cXYc=cca(Xc, Yc, simple=true)
@test cXYc.F[1]'*Cxc*cXYc.F[1]≈I
@test cXYc.F[2]'*Cyc*cXYc.F[2]≈I
Dc=cXYc.F[1]'*Cxyc*cXYc.F[2]
@test norm(Dc-Diagonal(Dc))+1. ≈ 1.
@test cXYc==cCc

# Method (3) real
# canonical correlation analysis of the average covariance and cross-covariance
k=10
Xset=[genDataMatrix(n, t) for i=1:k]
Yset=[genDataMatrix(n, t) for i=1:k]

c=cca(Xset, Yset)

# ... selecting subspace dimension allowing an explained variance = 0.9
c=cca(Xset, Yset; eVar=0.9)

# ... subtracting the mean from the matrices in Xset and Yset
c=cca(Xset, Yset; meanX=nothing, meanY=nothing, eVar=0.9)

# cca on the average of the covariance and cross-covariance matrices
# computed along dims 1
c=cca(Xset, Yset; dims=1, eVar=0.9)

# name of the filter
c.name

using Plots
# plot regularized accumulated eigenvalues
plot(c.arev)

# plot the original covariance and cross-covariance matrices
# and their transformed counterpart
 CxyMax=maximum(abs.(Cxy));
 h1 = heatmap(Cxy, clim=(-CxyMax, CxyMax), title="Cxy", yflip=true, c=:bluesreds);
 D=cC.F[1]'*Cxy*cC.F[2];
 Dmax=maximum(abs.(D));
 h2 = heatmap(D, clim=(0, Dmax), title="F1'CxyF2", yflip=true, c=:amp);
 CxMax=maximum(abs.(Cx));
 h3 = heatmap(Cx, clim=(-CxMax, CxMax), title="Cx", yflip=true, c=:bluesreds);
 h4 = heatmap(cC.F[1]'*Cx*cC.F[1], clim=(0, 1), title="F1'CxF1", yflip=true, c=:amp);
 CyMax=maximum(abs.(Cy));
 h5 = heatmap(Cy, clim=(-CyMax, CyMax), title="Cy", yflip=true, c=:bluesreds);
 h6 = heatmap(cC.F[2]'*Cy*cC.F[2], clim=(0, 1), title="F2'CyF2", yflip=true, c=:amp);
 📈=plot(h3, h5, h1, h4, h6, h2, size=(800,400))
# savefig(📈, homedir()*"\\Documents\\Code\\julia\\Diagonalizations\\docs\\src\\assets\\FigCCA.png")

# Method (3) complex
# canonical correlation analysis of the average covariance and cross-covariance
k=10
Xcset=[genDataMatrix(ComplexF64, n, t) for i=1:k]
Ycset=[genDataMatrix(ComplexF64, n, t) for i=1:k]

cc=cca(Xcset, Ycset)

# ... selecting subspace dimension allowing an explained variance = 0.9
cc=cca(Xcset, Ycset; eVar=0.9)

# ... subtracting the mean from the matrices in Xset and Yset
cc=cca(Xcset, Ycset; meanX=nothing, meanY=nothing, eVar=0.9)

# cca on the average of the covariance and cross-covariance matrices
# computed along dims 1
cc=cca(Xcset, Ycset; dims=1, eVar=0.9)

# name of the filter
cc.name

ext: jl


text:
using WNNMDenoise
using ImageTransformations, FileIO, ImageQualityIndexes, MAT
using Test

@testset "WNNMDenoise.jl" begin
    imgfile = joinpath(pkgdir(WNNMDenoise), "house.png")
    img = Float64.(load(imgfile) * 255)

    σₙ = 40
    noisy_img_file = joinpath(pkgdir(WNNMDenoise), "test", "house_AWGN_$(σₙ).mat")
    noisy_img = matopen(noisy_img_file) do io
        read(io, "N_Img")
    end

    @test assess_psnr(noisy_img, img, 255) ≈ 16.06 atol=1e-2

    img, noisy_img = restrict.((img, noisy_img))
    out = copy(noisy_img)
    f = WNNM(σₙ)
    f(out, noisy_img)

    println("PSNR: ", assess_psnr(out, img, 255))
    @test assess_psnr(out, img, 255) >= 28.18
end


text:
using SolidStateDetectors: getVe, getVh, Vl
using StaticArrays
using LinearAlgebra

config_dict = Dict(
    "High-purity germanium" => joinpath(@__DIR__,"../examples/example_config_files/ADLChargeDriftModel/drift_velocity_config.yaml"),
    "Silicon" => joinpath(@__DIR__,"../examples/example_config_files/ADLChargeDriftModel/drift_velocity_Si_300K_config.yaml")
)

geom_sigdigits(::Type{Int64})::Int = 12
geom_sigdigits(::Type{Float32})::Int = 6
geom_sigdigits(::Type{Float64})::Int = 12

for T in (Float32, Float64)
    @testset "Precision type: $(T)" begin

        for material in keys(config_dict)

            @testset "$material" begin
                config = config_dict[material]

                @testset "x-axis aligned with <100> axis" begin

                    cdm = ADLChargeDriftModel(config, T=T, phi110 = T(π/4)); # <100> aligned with x-axis
                    E = collect(T, 10 .^(0:0.1:6));

                    r100 = @SVector T[1,0,0]
                    r111 = @SVector T[sqrt(1/3), sqrt(1/3), sqrt(1/3)]

                    @testset "Electrons <100>" begin
                        ve100_cdm = broadcast(E -> abs(getVe(E * r100, cdm) ⋅ r100), E)
                        ve100_true = map(E -> Vl(E,cdm.electrons.axis100), E)
                        @test isapprox(ve100_cdm, ve100_true, rtol = 10.0^(-geom_sigdigits(T)))
                    end

                    @testset "Electrons <111>" begin
                        ve111_cdm = broadcast(E -> abs(getVe(E * r111, cdm) ⋅ r111), E)
                        ve111_true = map(E -> Vl(E,cdm.electrons.axis111), E)
                        @test isapprox(ve111_cdm, ve111_true, rtol = 10.0^(-geom_sigdigits(T)))
                    end

                    @testset "Holes <100>" begin
                        vh100_cdm = broadcast(E -> abs(getVh(E * r100, cdm) ⋅ r100), E)
                        vh100_true = map(E -> Vl(E,cdm.holes.axis100), E)
                        @test isapprox(vh100_cdm, vh100_true, rtol = 10.0^(-geom_sigdigits(T)))
                    end

                    @testset "Holes <111>" begin
                        vh111_cdm = broadcast(E -> abs(getVh(E * r111, cdm) ⋅ r111), E)
                        vh111_true = map(E -> Vl(E,cdm.holes.axis111), E)
                        @test isapprox(vh111_cdm, vh111_true, rtol = 10.0^(-geom_sigdigits(T)))
                    end
                end

                @testset "x-axis and <100> axis intersect at an angle of -30°" begin

                    cdm = ADLChargeDriftModel(config, T=T, phi110 = T(π/12)); # <110> and y-axis intersect at an angle of 15°
                    E = collect(T, 10 .^(0:0.1:6));

                    r100 = @SVector T[sqrt(3/4),-1/2,0]
                    r111 = @SVector T[(1 + sqrt(1/3))/2, (1 - sqrt(1/3))/2, sqrt(1/3)]

                    @testset "Electrons <100>" begin
                        ve100_cdm = broadcast(E -> abs(getVe(E * r100, cdm) ⋅ r100), E)
                        ve100_true = map(E -> Vl(E,cdm.electrons.axis100), E)
                        @test isapprox(ve100_cdm, ve100_true, rtol = 10.0^(-geom_sigdigits(T)))
                    end

                    @testset "Electrons <111>" begin
                        ve111_cdm = broadcast(E -> abs(getVe(E * r111, cdm) ⋅ r111), E)
                        ve111_true = map(E -> Vl(E,cdm.electrons.axis111), E)
                        @test isapprox(ve111_cdm, ve111_true, rtol = 10.0^(-geom_sigdigits(T)))
                    end

                    @testset "Holes <100>" begin
                        vh100_cdm = broadcast(E -> abs(getVh(E * r100, cdm) ⋅ r100), E)
                        vh100_true = map(E -> Vl(E,cdm.holes.axis100), E)
                        @test isapprox(vh100_cdm, vh100_true, rtol = 10.0^(-geom_sigdigits(T)))
                    end

                    @testset "Holes <111>" begin
                        vh111_cdm = broadcast(E -> abs(getVh(E * r111, cdm) ⋅ r111), E)
                        vh111_true = map(E -> Vl(E,cdm.holes.axis111), E)
                        @test isapprox(vh111_cdm, vh111_true, rtol = 10.0^(-geom_sigdigits(T)))
                    end
                end
            end
        end

        config = joinpath(@__DIR__,"../examples/example_config_files/ADLChargeDriftModel/drift_velocity_config_axes.yaml")

        @testset "All axes in one plane" begin

            cdm = ADLChargeDriftModel(config, T=T); # <110> and y-axis intersect at an angle of 15°
            E = collect(T, 10 .^(0:0.1:6));

            r100 = @SVector T[1,0,0]
            r111 = @SVector T[sqrt(1/3),sqrt(2/3),0]

            @testset "Electrons <100>" begin
                ve100_cdm = broadcast(E -> abs(getVe(E * r100, cdm) ⋅ r100), E)
                ve100_true = map(E -> Vl(E,cdm.electrons.axis100), E)
                @test isapprox(ve100_cdm, ve100_true, rtol = 10.0^(-geom_sigdigits(T)))
            end

            @testset "Electrons <111>" begin
                ve111_cdm = broadcast(E -> abs(getVe(E * r111, cdm) ⋅ r111), E)
                ve111_true = map(E -> Vl(E,cdm.electrons.axis111), E)
                @test isapprox(ve111_cdm, ve111_true, rtol = 10.0^(-geom_sigdigits(T)))
            end

            @testset "Holes <100>" begin
                vh100_cdm = broadcast(E -> abs(getVh(E * r100, cdm) ⋅ r100), E)
                vh100_true = map(E -> Vl(E,cdm.holes.axis100), E)
                @test isapprox(vh100_cdm, vh100_true, rtol = 10.0^(-geom_sigdigits(T)))
            end

            @testset "Holes <111>" begin
                vh111_cdm = broadcast(E -> abs(getVh(E * r111, cdm) ⋅ r111), E)
                vh111_true = map(E -> Vl(E,cdm.holes.axis111), E)
                @test isapprox(vh111_cdm, vh111_true, rtol = 10.0^(-geom_sigdigits(T)))
            end
        end
    end
end


text:
using Test, Astro

println("15.a Rise, Set, Transit")
# Meeus chapter 15, pages 103
# Venus 1988 March 20, at Boston (long 71.0833, lat 42.333)
# sidereal time at Greenwich in degrees should be 11h50m58.10 177d.74208

# coordinates of venus at 0h TD
# march 19 RA 40d.68021 Dec = 18d.04761
# march 20 RA 41.73129  Dec = 18d.44092
# march 21 RA 42.78204  Dec = 18d.82742

save_Lat = Astro.latitude
save_Long = Astro.longitude
set_longitude(deg2rad(71.0833))
set_latitude(deg2rad(42.3333))

println("using Longitude $(rad2deg(longitude))°" )
println("using Latitude $(rad2deg(latitude))°" )

println("Only low accuracy version currently available")

println("rise")
ut = cal_to_jd(1988, 3, 20)
ra = deg2rad(41.73129)
decl = deg2rad(18.44092)
# h0 is the 'standard' altitude, deg2rad(-0.5667) for stars/planets, deg2rad(-0.8333) for the sun
h0 = deg2rad(-0.5667)
jd_rise, jd_transit, jd_set = object_rts_low(ut, ra, decl, h0)

println("jd_rise $(fday_to_hms(jd_rise)), jd_transit $(fday_to_hms(jd_transit)), jd_set $(fday_to_hms(jd_set))")

#@test isapprox(jd_rise - ut, 0.51766, atol=0.1)

println("")

# restore lat and long
set_longitude(save_Long)
set_latitude(save_Lat)


text:
using GeometryTypes, Quaternions, GLVisualize, Reactive, GLAbstraction

function Quaternions.qrotation{T<:Real}(axis::Quaternion{T}, theta::T)
    ax = Vec{3,T}(axis.v1, axis.v2, axis.v3)
    qrotation(ax, theta)
end

function Quaternions.qrotation{T<:Real}(axis::Quaternion{T}, z::Quaternion{T}, theta::T)
    q = qrotation(axis, theta)
    q*z*conj(q)
end

function quaternion_belt_trick(timesignal)

    rad1 = 1f0 # outer circle
    rad2 = 0.2f0 # inner circle

    center1 = Point2f0(rad1,rad1) # left
    center2 = Point2f0(rad1+2f0,rad1) # right

    u1234 = Quaternion(0f0,0f0,1f0,0f0) # gets rotated
    axis  = Quaternion(0f0,0f0,0f0,1f0) # axis of qrotationdd
    mi  = 80
    mj  = 10
    rxs = zeros(Float32, mi+1,mj+1)
    rys = zeros(Float32, mi+1,mj+1)
    rzs = zeros(Float32, mi+1,mj+1)

    max_frames = 96
    dphi = 2*pi/(max_frames-1) # include both ends

    u1234_s = foldp(u1234, timesignal) do v0, _
        qrotation(axis, v0, Float32(dphi))
    end
    xyz = const_lift(plot_belts, rxs, rys, rzs, u1234_s, center1, center2, rad1, rad2)
    x = const_lift(getindex, xyz, 1)
    y = const_lift(getindex, xyz, 2)
    z = const_lift(getindex, xyz, 3)
    (x,y,z)
end

function plot_belts(rxs, rys, rzs, u1234, center1, center2, rad1, rad2)
    mi  = 80
    mj  = 10
    iis = 1f0:(mi+1f0)
    rs  = rad1-((rad1-rad2)*(iis-1f0)/mi)
    for ii=1:(mi+1)
        jjs   = 1f0:(mj+1f0)
        xs    = (rad2/mj)*((jjs-1f0)-(mj/2f0))
        ys    = sqrt(rs[ii]*rs[ii]-xs.*xs)
        theta = Float32(2*pi*(ii-1)/mi)

        for jj=1:(mj+1)
            efgh    = Quaternion(0f0, xs[jj], ys[jj], 0f0)
            wxyz    = qrotation(efgh, u1234, theta)
            rxs[ii,jj] = center1[1]+wxyz.v1
            rys[ii,jj] = center1[2]+wxyz.v2
            rzs[ii,jj] = wxyz.v3
        end
    end
    rxs, rys, rzs
end

if !isdefined(:runtests)
    window = glscreen()
    timesignal = loop(linspace(0f0, 1f0, 360))
end

xyz = quaternion_belt_trick(timesignal)
_view(visualize(xyz, :surface), window)

if !isdefined(:runtests)
    renderloop(window)
end


text:
include("../examples/lorenz63.jl")
using Test
using DiffEqSensitivity, OrdinaryDiffEq, ForwardDiff
using Zygote
function time_average_obj(s, n, u0)
	u = copy(u0)
	J_avg = 0.
	for i = 1:n
		u = lorenz63_ad(u, s, 1)
		J_avg += u[end]/n
	end
	return J_avg
end
	s = [10., 28., 8/3]
	m = 1
    n = 200
    n_runup = 5000
	n_samples = 1000
    u0 = rand(3)
	u0 = lorenz63_ad(u0, s, n_runup)
	dJds = zeros(3,n_samples)
	for i = 1:n_samples
		u0 .= lorenz63_ad(u0, s, 1)
		dJds[:,i] .= ForwardDiff.gradient(s->time_average_obj(
										s,n,u0), s)
	end
	@show sum(dJds, dims=2)/n_samples

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
export DepthwiseConvDims

"""
    DepthwiseConvDims

Concrete subclass of `ConvDims` for a depthwise convolution.  Differs primarily due to
characterization by C_in, C_mult, rather than C_in, C_out.  Useful to be separate from
DenseConvDims primarily for channel calculation differences.
"""
struct DepthwiseConvDims{N,S,P,D,F} <: ConvDims{N,S,P,D,F}
    I::NTuple{N, Int}
    K::NTuple{N, Int}
    C_in::Int
    C_mult::Int
end

# Getters for the fields
input_size(c::DepthwiseConvDims) = c.I
kernel_size(c::DepthwiseConvDims) = c.K
channels_in(c::DepthwiseConvDims) = c.C_in
channels_out(c::DepthwiseConvDims) = c.C_in * channel_multiplier(c)
channel_multiplier(c::DepthwiseConvDims) = c.C_mult

# Convenience wrapper to create DepthwiseConvDims objects
function DepthwiseConvDims(x_size::NTuple{M}, w_size::NTuple{M};
                           stride=1, padding=0, dilation=1, flipkernel::Bool=false) where M
    # Do common parameter validation
    stride, padding, dilation = check_spdf(x_size, w_size, stride, padding, dilation)

    # Ensure channels are equal
    if x_size[end-1] != w_size[end]
        xs = x_size[end-1]
        ws = w_size[end]
        throw(DimensionMismatch("Input channels must match! ($xs vs. $ws)"))
    end

    return DepthwiseConvDims{
        M - 2,
        stride,
        padding,
        dilation,
        flipkernel
    }(
        # Image spatial size
        x_size[1:end-2],

        # Kernel spatial size
        w_size[1:end-2],

        # Input channels
        x_size[end-1],

        # Channel multiplier
        w_size[end-1],
    )
end

# Auto-extract sizes and just pass those directly in
function DepthwiseConvDims(x::AbstractArray, w::AbstractArray; kwargs...)
    if ndims(x) != ndims(w)
        throw(DimensionMismatch("Rank of x and w must match! ($(ndims(x)) vs. $(ndims(w)))"))
    end
    return DepthwiseConvDims(size(x), size(w); kwargs...)
end

# Useful for constructing a new DepthwiseConvDims that has only a few elements different
# from the original progenitor object.
function DepthwiseConvDims(c::DepthwiseConvDims; N=spatial_dims(c), I=input_size(c), K=kernel_size(c),
                           C_in=channels_in(c), C_m=channel_multiplier(c), S=stride(c),
                           P=padding(c), D=dilation(c), F=flipkernel(c))
    return DepthwiseConvDims{N, S, P, D, F}(I, K, C_in, C_m)
end

# This one is basically the same as for DenseConvDims, we only change a few lines for kernel channel count
function check_dims(x::NTuple{M}, w::NTuple{M}, y::NTuple{M}, cdims::DepthwiseConvDims) where {M}
    # First, check that channel counts are all correct:
    @assert x[end-1] == channels_in(cdims) DimensionMismatch("Data input channel count ($(x[end-1]) vs. $(channels_in(cdims)))")
    @assert y[end-1] == channels_out(cdims) DimensionMismatch("Data output channel count ($(y[end-1]) vs. $(channels_out(cdims)))")
    @assert w[end-1] == channel_multiplier(cdims) DimensionMismatch("Kernel multiplier channel count ($(w[end-1]) vs. $(channel_multiplier(cdims))")
    @assert w[end] == channels_in(cdims) DimensionMismatch("Kernel input channel count ($(w[end]) vs. $(channels_in(cdims)))")

    # Next, check that the spatial dimensions match up
    @assert x[1:end-2] == input_size(cdims) DimensionMismatch("Data input spatial size ($(x[1:end-2]) vs. $(input_size(cdims)))")
    @assert y[1:end-2] == output_size(cdims) DimensionMismatch("Data output spatial size ($(y[1:end-2]) vs. $(output_size(cdims)))")
    @assert w[1:end-2] == kernel_size(cdims) DimensionMismatch("Kernel spatial size ($(w[1:end-2]) vs. $(kernel_size(cdims)))")

    # Finally, check that the batch size matches
    @assert x[end] == y[end] DimensionMismatch("Batch size ($(x[end]) vs. $(y[end]))")
end


text:
using Plots
using LLLplus
using Printf
using Random
using LinearAlgebra
using DoubleFloats
using Quadmath
gr();

#include("../src/l2.jl")
include("lrtest.jl")

# getIntType is used to indicate what type we want the unimodular integer
# matrix to have. We need to add methods of getIntType to handle the
# DoubleFloats and Quadmath, since they are external packages that LLLplus
# doesn't know about. If this is type piracy, I think it's acceptable.
import LLLplus.getIntType
getIntType(Td::Type{Tr}) where {Tr<:Float128} = Int128
getIntType(Td::Type{Tr}) where {Tr<:Double64} = Int128

lrtest(40,2 .^[3],[100],[Int32,Int64,Int128,Float32,Float64,Double64,Float128,BigInt,BigFloat],"rand")
savefig("perfVsDataType.png")

# lrtest(40,2 .^[1:8;],[1],[Float64],"randn")
# savefig("perfVsNfloat64.png")


text:
using FastRunningMedian, Test, DataStructures, JLD2
import Statistics

"""
    check_health(mf::MedianFilter)

Check that all pointers point at a thing that again points back at them. Also check that low_heap is smaller than high_heap.

Debug Function for MedianFilter.
"""
function check_health(mf::MedianFilter)
    for k in 1:length(mf.heap_pos)
        current_heap, current_heap_ind = mf.heap_pos[k]
        if current_heap == true
            a = mf.low_heap[current_heap_ind][2] - mf.heap_pos_offset
            # println(k, " ?= ", a)
            @assert k == a
        else
            a = mf.high_heap[current_heap_ind][2] - mf.heap_pos_offset
            # println(k, " ?= ", a)
            @assert k == a
        end
    end

    if length(mf) >= 2
        @assert first(mf.low_heap) <= first(mf.high_heap)
    end
end

println("running tests...")

@testset "FastRunningMedian Tests" begin

    @testset "Stateful API Tests" begin

        @testset "Grow and Shrink Fuzz" begin
            # load test cases
            @load "fixtures/grow_shrink.jld2" grow_shrink_fixtures

            function grow_and_shrink_test(values, expected_medians)
                N = length(values)
                mf = MedianFilter(values[1], N)
                check_health(mf)
                @assert median(mf) == expected_medians[1]
                #grow phase
                for i in 2:N
                    grow!(mf, values[i])
                    check_health(mf)
                    @assert median(mf) == expected_medians[i]
                end
                # shrink phase
                for i in N+1:2N-1
                    shrink!(mf)
                    check_health(mf)
                    @assert median(mf) == expected_medians[i]
                end
                @assert length(mf) == 1
            end

            for fixture in grow_shrink_fixtures
                grow_and_shrink_test(fixture...)
            end
        end

        @testset "Roll Fuzz" begin
            #load test cases
            @load "fixtures/roll.jld2" roll_fixtures

            function roll_test(initial_values, roll_values, expected_medians)
                window_size = length(initial_values)
                mf = MedianFilter(initial_values[1], window_size)
                for i in 2:window_size
                    grow!(mf, initial_values[i])
                end
                @assert length(mf) == window_size
                for i in 1:length(roll_values)
                    mf_median = roll!(mf, roll_values[i])
                    check_health(mf)
                    @assert mf_median == expected_medians[i]
                end
            end

            for fixture in roll_fixtures
                roll_test(fixture...)
            end
        end

        @testset "Roll does work with window size 1" begin
            mf = MedianFilter(1., 1)
            @test 2. == roll!(mf, 2.)
            @test 1. == roll!(mf, 1.)
        end

        @testset "Grow! does not grow beyond capacity" begin
            mf = MedianFilter(1., 3)
            grow!(mf, 2.)
            grow!(mf, 3.)
            @test_throws ErrorException grow!(mf, 4.)
        end

        @testset "shrink! below 1-element errors" begin
            mf = MedianFilter(1., 4)
            @test_throws ErrorException shrink!(mf)
        end

        @testset "can only roll when capacity is exactly met" begin
            mf = MedianFilter(1., 3)
            grow!(mf, 2.)
            @test_throws ErrorException roll!(mf, 3.)
            grow!(mf, 3.)
            roll!(mf, 4.)
            shrink!(mf)
            @test_throws ErrorException roll!(mf, 5.)
        end
    end

    @testset "High Level API Tests" begin
        # Desired API
        # running_median(input::Array{T, 1}, window_size::Integer, tapering=:sym) where T <: Real
        # taperings:
        # :symmetric or :sym (window symmetric around returned point, length N-1 if even window, N if odd)
        # :asymmetric or :asym (window full length to one side, length length N+W-1 if odd W, N-1+W if even window)
        # :asymmetric_truncated or :asymtrunc (same as asymmetric, but truncated at ends to size of symmetric)
        # :none or :no (only full length window used, length N-W+1)
        #
        # all these taperings are symmetrical in that they behave the same at each end of the array, only mirrored

        @testset "Basic API examples" begin
            @test_throws ErrorException running_median(zeros(0), 1)
            @test running_median([1.], 1) == [1.]
            @test running_median([1., 2., 3.], 1) == [1., 2., 3.]
            @test running_median([1., 4., 2., 1.], 3) == [1., 2., 2., 1.]
            @test running_median([1, 4, 2, 1], 3) == [1, 2, 2, 1]
            @test running_median([1, 4, 2, 1], 3, :asym) == [1, 2.5, 2, 2, 1.5, 1]
            @test running_median([1., 4., 2., 1.], 3, :sym) == [1., 2., 2., 1.]
            @test running_median([1., 4., 2., 1.], 3, :symmetric) == [1., 2., 2., 1.]
            @test running_median([1., 4., 2., 1.], 3, :asym) == [1., 2.5, 2., 2., 1.5, 1.]
            @test running_median([1., 4., 2., 1.], 3, :asymmetric) == [1., 2.5, 2., 2., 1.5, 1.]
            @test running_median([1., 4., 2., 1.], 3, :asym_trunc) == [2.5, 2., 2., 1.5]
            @test running_median([1., 4., 2., 1.], 3, :asymmetric_truncated) == [2.5, 2., 2., 1.5]
            @test running_median([1., 4., 2., 1.], 3, :no) == [2., 2.]
            @test running_median([1., 4., 2., 1.], 3, :none) == [2., 2.]
            @test running_median([1., 2., 1., 2., 1., 3.], 101) == [1., 1., 1., 2., 2., 3.]
            @test running_median([1., 1., 2., 1., 1., 1., 1., 1., 2., 1.], 99) ==
                [1., 1., 1., 1., 1., 1., 1., 1., 1., 1.]
        end

        @testset "Compare to Naive Symmetric Median" begin
            @load "fixtures/symmetric.jld2" fixtures
            for fixture in fixtures
                @test fixture[3] == running_median(fixture[1], fixture[2], :sym)
            end
        end

        @testset "Compare to Naive Asymmetric Median" begin
            @testset "Float Input" begin
                @load "fixtures/asymmetric.jld2" fixtures
                for fixture in fixtures
                    @test fixture[3] == running_median(fixture[1], fixture[2], :asym)
                end
            end

            @testset "Int Input" begin
                @load "fixtures/asymmetric_int.jld2" fixtures
                for fixture in fixtures
                    @test fixture[3] == running_median(fixture[1], fixture[2], :asym)
                end
            end
        end

        @testset "Compare to Untapered Median from RollingFunctions" begin
            @load "fixtures/untapered.jld2" untapered_fixtures
            for fixture in untapered_fixtures
                @test fixture[3] == running_median(fixture[1], fixture[2], :none)
            end
        end

        @testset "Compare to Naive Asymmetric Truncated Median" begin
            @load "fixtures/asym_trunc.jld2" asym_trunc_fixtures
            for fixture in asym_trunc_fixtures
                @test fixture[3] == running_median(fixture[1], fixture[2], :asym_trunc)
            end
        end

    end
end # all tests

max_stars_repo_name: Firionus/FastRunningMedian

max_issues_repo_name: Firionus/FastRunningMedian


text:
using Test
using LinearAlgebra
using Nystrom

@testset "Area/volume" begin
    @testset "Cube" begin
        # generate a mesh
        (lx,ly,lz) = widths = (1.,1.,2.)
        Ω,M = @gmsh begin
            Geometry.clear_entities!()
            gmsh.model.occ.addBox(0,0,0,lx,ly,lz)
            gmsh.model.occ.synchronize()
            gmsh.model.mesh.generate(3)
            Ω = GmshSDK.domain(dim=3)
            M = GmshSDK.meshgen(Ω,dim=3)
            return Ω,M
        end
        ∂Ω = boundary(Ω)
        mesh  = NystromMesh(view(M,∂Ω),order=1)
        A     = 2*(lx*ly + lx*lz + ly*lz)
        @test A ≈ sum(qweights(mesh))
        # generate a Nystrom mesh for volume
        mesh  = NystromMesh(view(M,Ω),order=1)
        V     = prod(widths)
        # sum only weights corresponding to tetras
        @test V ≈ sum(qweights(mesh))
    end
    @testset "Sphere" begin
        # create a mesh using GmshSDK package
        r = 0.5
        Geometry.clear_entities!()
        Ω,M   = @gmsh begin
            GmshSDK.set_meshsize(0.1)
            Geometry.clear_entities!()
            gmsh.model.occ.addSphere(0,0,0,r)
            gmsh.model.occ.synchronize()
            gmsh.model.mesh.generate(3)
            Ω = GmshSDK.domain(dim=3)
            M = GmshSDK.meshgen(Ω,dim=3)
            return Ω,M
        end
        Γ     = boundary(Ω)
        nmesh = NystromMesh(M,Γ,order=4) # NystromMesh of surface Γ
        area  = sum(qweights(nmesh))
        @test isapprox(area,4*π*r^2,atol=5e-2)
        nmesh = NystromMesh(M,Ω,order=4) # Nystrom mesh of volume Ω
        volume  = sum(qweights(nmesh))
        @test isapprox(volume,4/3*π*r^3,atol=1e-2)
        # create a mesh using ParametricSurfaces package
        Geometry.clear_entities!()
        ent   = ParametricSurfaces.Sphere(radius=r)
        Ω     = Geometry.Domain(ent)
        Γ     = boundary(Ω)
        M     = ParametricSurfaces.meshgen(Γ,(4,4))
        nmesh = NystromMesh(view(M,Γ),order=5)
        area   = sum(qweights(nmesh))
        @test isapprox(area,4*π*r^2,atol=1e-5)
    end
    @testset "Circle" begin
        r = rx = ry = 0.5
        Geometry.clear_entities!()
        Ω,M = @gmsh begin
            Geometry.clear_entities!()
            gmsh.model.occ.addDisk(0,0,0,rx,ry)
            gmsh.model.occ.synchronize()
            gmsh.model.mesh.generate(2)
            Ω = GmshSDK.domain(dim=2)
            M = GmshSDK.meshgen(Ω,dim=2)
            return Ω,M
        end
        Γ    = boundary(Ω)
        mesh = NystromMesh(view(M,Ω),order=2)
        A    = π*r^2
        # test area
        @test isapprox(A,sum(qweights(mesh));atol=1e-2)
        # test perimeter
        mesh = NystromMesh(view(M,Γ),order=2)
        P    = 2π*r
        @test isapprox(P,sum(qweights(mesh));atol=1e-2)
        Geometry.clear_entities!()
        geo   = ParametricSurfaces.Circle(radius=r)
        Ω     = Geometry.Domain(geo)
        Γ     = boundary(Ω)
        M     = ParametricSurfaces.meshgen(Γ,(100,))
        nmesh = NystromMesh(M,Γ,order=5)
        @test isapprox(sum(qweights(nmesh)),π,atol=1e-10)
    end
end

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom


text:
module SourceTest

using HPAT
#HPAT.DomainPass.set_debug_level(3)
#HPAT.DataTablePass.set_debug_level(3)

@acc hpat function stest(file_name)
    A = DataSource(Vector{Float64},HDF5,"/A", file_name)
    return sum(A)
end

@acc hpat function stest2(file_name)
    A = DataSource(Matrix{Float64},HDF5,"/A2", file_name)
    return sum(A)
end

using HDF5
using MPI
using Base.Test

function main()

    if MPI.Comm_rank(MPI.COMM_WORLD)==0
        h5write("A.hdf5", "/A", [1.0,2.0,3.0,4.0])
        h5write("A.hdf5", "/A2", [1.0 2.0 3.0; 4.0 5.0 6.0])
    end

    @test_approx_eq SourceTest.stest("A.hdf5") 10.0
    @test_approx_eq SourceTest.stest2("A.hdf5") 21.0
    if MPI.Comm_rank(MPI.COMM_WORLD)==0 rm("A.hdf5") end
end

end

SourceTest.main()


text:
# This file is a part of ROOTFramework.jl, licensed under the MIT License (MIT).

using ROOTFramework

using Compat.Test

import JSON

@testset "TTree I/O" begin
    hist = TH1D()
    JSON.parse(rootjson(hist))
end


text:
using DigitalComm
using PGFPlotsX

function monteCarlo(snr,mcs,nbSymb)
	# Number of bits
	nbBits	  = nbSymb * Int(log2(mcs));
	# --- Binary sequence generation
	bitSeq	  = genBitSequence(nbBits);
	# --- QPSK mapping
	qamSeq	  = bitMappingQAM(mcs,bitSeq);
	# ----------------------------------------------------
	# --- Channel
	# ----------------------------------------------------
	#  --- AWGN
	# Theoretical power is 1 (normalized constellation)
	qamNoise,  = addNoise(qamSeq,snr,1);
	# ----------------------------------------------------
	# --- Rx Stage: SRRC
	# ----------------------------------------------------
	# --- Binary demapper
	bitDec	= bitDemappingQAM(mcs,qamNoise);
	# --- Error counter
	nbE		= sum(xor.(bitDec,bitSeq));
	# --- Return Error and bits
	return (nbE,nbBits);
end

function doPlot(snrVect,ber,qamVect)
	a = 0;
	@pgf a = Axis({
				   ymode	  = "log",
				   height      ="3in",
				   width       ="4in",
				   grid,
				   xlabel      = "SNR [dB]",
				   ylabel      = "Bit Error Rate ",
				   ymax 	   = 1,
				   ymin 	   = 10.0^(-5),
				   title       = "AWGN BER for QAM",
				   legend_style="{at={(0,0)},anchor=south west,legend cell align=left,align=left,draw=white!15!black}"
				   },
				  Plot({color="red",mark="square*"},Table([snrVect,ber[1,:]])),
				  LegendEntry("QPSK"),
				  Plot({color="green",mark="*"},Table([snrVect,ber[2,:]])),
				  LegendEntry("16-QAM"),

				  Plot({color="purple",mark="triangle*"},Table([snrVect,ber[3,:]])),
				  LegendEntry("64-QAM"),
				  Plot({color="blue",mark="diamond*"},Table([snrVect,ber[4,:]])),
				  LegendEntry("256-QAM"),
				  );
	# ---  Adding theoretical curve
	snrLin  = (10.0).^(snrVect/10)
	for qamScheme = qamVect
		ebNo 	= snrLin / log2(qamScheme);
		# This approximation is only valid for high SNR (one symbol error is converted to one bit error with Gray coding).
		berTheo	  = 4 * ( 1 - 1 / sqrt(qamScheme)) / log2(qamScheme) * qFunc.(sqrt.( 2*ebNo * 3 * log2(qamScheme) / (2*(qamScheme-1)  )));
		@pgf push!(a,Plot({color="black"},Table([snrVect,berTheo])));
	end
	display(a);
end

function main()
	# --- Parameters
	nbIt		    = 10000;			  # Number of iterations
	nbSymb	        = 1024;			  # Number of symbols per iterations
	mcs		        = [4,16,64,256];  # Constellation size
	snrRange	    = (-1:26);		  # SNR, expressed in dB
	# --- Init performance metrics
	nbSNR			= length(snrRange);
	ber			    = zeros(Float64,length(mcs),nbSNR);
	for iMcs = 1 : 1  : length(mcs)
		for iSNR = 1 : 1 : nbSNR
			# --- Create BER counters
			nbE	  = 0;
			nbB	  = 0;
			for iN = 1 : 1 : nbIt
				# --- Elementary MC call
				# Corresponds to a given SNR and a given iteration
				# As we are ergodic in AWGN, it is only nbSymb*nbIt that matters for BER computation
				(a,b) = monteCarlo(snrRange[iSNR],mcs[iMcs],nbSymb);
				# --- Update counters
				nbE	+= a;	  # Increment errors
				nbB	+= b;	  # Increment bit counters
			end
			ber[iMcs,iSNR] = nbE / nbB;
		end
	end
	# --- Plotting routine
	doPlot(snrRange,ber,mcs);
end


text:
import java.security.KeyPairGenerator;
import java.security.KeyPair;
import java.security.PublicKey;
import java.security.PrivateKey;
import javax.crypto.Cipher;

public class CryptoLibrary {
    public static byte[] encrypt(PublicKey key, byte[] plaintext) throws Exception {
	return new byte[plaintext.length];
    }

    public static byte[] decrypt(PrivateKey key, byte[] ciphertext) throws Exception {
	return new byte[ciphertext.length];
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
"""
Smooth inverse frequency (SIF) structure for document embedding
using word vectors.

# References
  * [Arora et al. ICLR 2017, "A simple but tough-to-beat baseline for sentence embeddings"]
    (https://openreview.net/pdf?id=SyK00v5xx)
"""
struct SIFEmbedder{T,S} <: WordVectorsEmbedder{T,S}
    embeddings::EmbeddingsLibrary{T,S}
    lexicon::OrderedDict{S, Int}
    alpha::T
    config::NamedTuple
end

SIFEmbedder(embeddings::EmbeddingsLibrary{T,S},
            config;
            lexicon=OrderedDict{S, Int}(),
            alpha=DEFAULT_SIF_ALPHA,
            kwargs...
            ) where {T,S} =
    SIFEmbedder(embeddings, lexicon, T(alpha), config)

# Dimensionality function
function dimensionality(embedder::SIFEmbedder)
    return size(embedder.embeddings)[1]
end

# Sentence embedding function
function sentences2vec(embedder::SIFEmbedder,
                       document_embedding::Vector{Matrix{T}};
                       embedded_words::Vector{Vector{S}}=[String[]],
                       kwargs...) where {T,S}
    isempty(document_embedding) && return zeros(T, dimensionality(embedder), 0)
    sif(document_embedding,
        embedder.lexicon,
        embedded_words,
        embedder.alpha,
        dimensionality(embedder))
end

"""
    sif(document_embedding, lexicon, embedded_words, alpha, dim)

Implementation of sentence embedding principled on subtracting the paragraph vector i.e.
principal vector of a sentence from the sentence's word embeddings.
"""
#TODO(Corneliu): Make the calculation of `alpha` automatic using some heuristic
function sif(document_embedding::Vector{Matrix{T}},
             lexicon::OrderedDict{String, Int},
             embedded_words::Vector{Vector{S}},
             alpha::T,
             dim::Int
            ) where {T<:AbstractFloat, S<:AbstractString}
    L = sum(values(lexicon))
    n = length(document_embedding)  # number of sentences in document
    X = zeros(T, dim, n)  # new document embedding
    # Loop over sentences
    for (i, s) in enumerate(document_embedding)
        p = [get(lexicon, word, eps(T))/L for word in embedded_words[i]]
        W = size(s,2)  # no. of words
        @inbounds for w in 1:W
            X[:,i] += 1/(length(s)) * (alpha/(alpha+p[w]) .* s[:,w])
        end
    end
    local u::Vector{T}
    try
        u₀, _, _ = tsvd(X, 1)
        u = vec(u₀)
    catch
        u₀, _, _ = svd(X)
        u =u₀[:, 1]
    end
    @inbounds @simd for i in 1:n
        X[:,i] -= (u*u') * X[:,i]
    end
    return X
end


text:
# Unit test for EOP C04 file reading
@testset "EOP C04 file reading" begin
	data = loadeop(joinpath(dirname(@__DIR__),"test","input","eop_data.c04"));
	@test size(data) == (11,14)
	@test data[:datetime][1] == DateTime(1962,1,1)
	@test data[:datetime][end] == DateTime(1962,1,11)
	@test sum(data[:x]) ≈ -0.297387
	@test data[:dYErr][end] ≈ 0.002000
end


text:
using Documenter, Juniper

makedocs(
    # See https://github.com/JuliaDocs/Documenter.jl/issues/868
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true", mathengine = Documenter.MathJax()),
    # See https://github.com/JuliaOpt/JuMP.jl/issues/1576
    strict = true,
    sitename = "Juniper",
    pages = [
        "Home" => "index.md",
        "Options" => "options.md",
        "Extras" => "extras.md",
#        "Developer" => [],
#        "Library" => "library.md"
    ]
)

deploydocs(
    repo = "github.com/lanl-ansi/Juniper.jl.git",
)


text:
"""
    (dvec_ngwp, BS_ngwp) = ngwp_bestbasis(dmatrix::Array{Float64,3}, GP_star::GraphPart;
                                          cfspec::Any = 1.0, flatten::Any = 1.0,
                                          j_start::Int = 1, j_end::Int = size(dmatrix, 2),
                                          useParent::Bool = true)

Select the best basis from the matrix of NGWP expansion coefficients.

# Input Arguments
- `dmatrix::Array{Float64,3}`: the matrix of expansion coefficients
- `GP_star::GraphPart`: an input GraphPart object of the dual graph
- `cfspec::Any`: the specification of cost functional to be used (default = 1.0,
    i.e., 1-norm)
- `flatten::Any`: the method for flattening vector-valued data to scalar-valued
    data (default = 1.0, i.e, 1-norm)
- `useParent::Bool`: the flag to indicate if we update the selected best basis
    subspace to the parent when parent and child have the same cost (default = false)

#  Output Arguments
- `dvec_ngwp::Matrix{Float64}`: the vector of expansion coefficients corresponding
    to the NGWP best basis
- `BS_ngwp::BasisSpec`: a BasisSpec object which specifies the NGWP best basis
"""
function ngwp_bestbasis(dmatrix::Array{Float64,3}, GP_star::GraphPart;
                            cfspec::Any = 1.0, flatten::Any = 1.0,
                            j_start::Int = 1, j_end::Int = size(dmatrix, 2),
                            useParent::Bool = false)
    dvec_ngwp, BS_ngwp = ghwt_c2f_bestbasis(dmatrix, GP_star; cfspec = cfspec,
                            flatten = flatten, j_start = j_start, j_end = j_end,
                            useParent = useParent)
    BS_ngwp.description = "NGWP Best Basis"
    return dvec_ngwp, BS_ngwp
end


text:
function numcollatz(n)
	numiters = 0
	while n > 1
		#print("$n ")
		if n % 2 == 0
			n = div(n,2)
		else
			n = 3n + 1
		end
		#println("$n")
		numiters += 1
	end
	return numiters
end

function main()
	N = 1_000_000

	bests = 0
	bestl = 0
	i = 0
	while i <= N
		i += 1
		l = numcollatz(i)
		if l > bestl
			println("$i $l")
			bestl = l
			bests = i
		end
	end
	println("start=$bests length=$bestl")
end

main()

max_stars_repo_name: cmey/euler

max_issues_repo_name: cmey/euler

max_forks_repo_name: cmey/euler


text:
export SingPoint,
    LogSingPoint,
    AlgebraicSingPoint,
    SingularDiagonal

@deprecate DiagonallySingular SingularDiagonal

"""
A `Property` type describes a property of the integrand that is independent
from the other aspects of the integral (domain, measure).

The main examples are singularities of the integrand function, i.e., smoothness
properties.
"""
abstract type Property end

"No property of the integrand is known."
struct NoProperty <: Property end

"The supertype of all kinds of singularities of an integrand."
abstract type Singularity <: Property end

"No singularity of the integrand is known."
struct NoSingularity <: Singularity end

"Supertype of all kinds of point singularities."
abstract type PointSingularity <: Singularity end

point(s::PointSingularity) = s.point

"An unspecified singularity occurring at a specific point."
struct SingPoint{T} <: PointSingularity
    point   ::  T
end

"An unspecified singularity occurring at a corner point."
struct SingCorner{T} <: PointSingularity
    point   ::  T
end

"A logarithmic singularity at a specific point."
struct LogSingPoint{T} <: PointSingularity
    point   ::  T
end

"An algebraic singularity at a specific point."
struct AlgebraicSingPoint{O,T} <: PointSingularity
    point   ::  T
    order   ::  O
end

order(s::AlgebraicSingPoint) = s.order

"Supertype of singularities along a curve."
abstract type CurveSingularity <: Singularity end

"For 2D integrands, a singularity along the line `x=y`."
struct SingularDiagonal <: CurveSingularity
end

"""
A list of breakpoints of the integrand, for example because the integrand
is piecewise smooth.
"""
struct Waypoints{P} <: Property
    points  ::  P
end

points(P::Waypoints) = P.points


text:
include("test_kernel.jl")
include("test_poly.jl")
include("test_basis.jl")
include("test_basis_set.jl")
include("test_sampling.jl")
include("test_augmentation.jl")
include("test_composite.jl")
#include("test_spr.jl") # Fails


text:
function _container_spec(V::DataType, ax...; kwargs...)

    parameters = get(kwargs, :parameters, true)

    # While JuMP fixes the isassigned problems
     # While JuMP fixes the isassigned problems
    #=
    if parameters
            cont = JuMP.Containers.DenseAxisArray{PGAE{V}}(undef, ax...)
            _remove_undef!(cont.data)
        return cont
    else
            cont = JuMP.Containers.DenseAxisArray{GAE{V}}(undef, ax...)
            _remove_undef!(cont.data)
        return cont
    end
    =#

    if parameters
        return JuMP.Containers.DenseAxisArray{PGAE{V}}(undef, ax...)
    else
        return JuMP.Containers.DenseAxisArray{GAE{V}}(undef, ax...)
    end

end

function _canonical_model_init(bus_numbers::Vector{Int64},
                              optimizer::Union{Nothing,JuMP.OptimizerFactory},
                              transmission::Type{S},
                              time_range::UnitRange{Int64}; kwargs...) where {S <: PM.AbstractPowerFormulation}

    parameters = get(kwargs, :parameters, true)

    jump_model = _pass_abstract_jump(optimizer; kwargs...)
    V = JuMP.variable_type(jump_model)

    ps_model = CanonicalModel(jump_model,
                            Dict{Symbol, JuMP.Containers.DenseAxisArray}(),
                            Dict{Symbol, JuMP.Containers.DenseAxisArray}(),
                            zero(JuMP.GenericAffExpr{Float64, V}),
                            Dict{Symbol, JuMP.Containers.DenseAxisArray}(:var_active => _container_spec(V, bus_numbers, time_range; kwargs...),
                                                                         :var_reactive => _container_spec(V, bus_numbers, time_range; kwargs...)),
                            parameters ? Dict{Symbol,JuMP.Containers.DenseAxisArray}() : nothing,
                            Dict{Symbol,Array{InitialCondition}}(),
                            nothing);

    return ps_model

end

function _canonical_model_init(bus_numbers::Vector{Int64},
                               optimizer::Union{Nothing,JuMP.OptimizerFactory},
                               transmission::Type{S},
                               time_range::UnitRange{Int64}; kwargs...) where {S <: PM.AbstractActivePowerFormulation}

    parameters = get(kwargs, :parameters, true)

    jump_model = _pass_abstract_jump(optimizer; kwargs...)
    V = JuMP.variable_type(jump_model)

    ps_model = CanonicalModel(jump_model,
                              Dict{Symbol, JuMP.Containers.DenseAxisArray}(),
                              Dict{Symbol, JuMP.Containers.DenseAxisArray}(),
                              zero(JuMP.GenericAffExpr{Float64, V}),
                              Dict{Symbol, JuMP.Containers.DenseAxisArray}(:var_active => _container_spec(V, bus_numbers, time_range; kwargs...)),
                              parameters ? Dict{Symbol,JuMP.Containers.DenseAxisArray}() : nothing,
                              Dict{Symbol,Array{InitialCondition}}(),
                              nothing);

    return ps_model

end

function  build_canonical_model(transmission::Type{T},
                                devices::Dict{Symbol, DeviceModel},
                                branches::Dict{Symbol, DeviceModel},
                                services::Dict{Symbol, ServiceModel},
                                system::PSY.System,
                                optimizer::Union{Nothing,JuMP.OptimizerFactory}=nothing;
                                kwargs...) where {T <: PM.AbstractPowerFormulation}

time_range = 1:system.time_periods
bus_numbers = [b.number for b in system.buses]

ps_model = _canonical_model_init(bus_numbers, optimizer, transmission, time_range; kwargs...)

# Build Injection devices
for mod in devices
    construct_device!(ps_model, mod[2], transmission, system, time_range; kwargs...)
end

# Build Network
construct_network!(ps_model, transmission, system, time_range; kwargs...)

# Build Branches
for mod in branches
    construct_device!(ps_model, mod[2], transmission, system, time_range; kwargs...)
end

#Build Service
for mod in services
    construct_service!(ps_model, mod[2], transmission, system, time_range; kwargs...)
end

# Objective Function
JuMP.@objective(ps_model.JuMPmodel, Min, ps_model.cost_function)

return ps_model

end


text:
using Test
import AstroBase

using AstroTime: Epoch, TDBEpoch, SECONDS_PER_DAY, j2000, seconds, value, julian_twopart
using AstroBase.Ephemerides
using AstroBase.Bodies
using AstroBase.Constants: astronomical_unit
using ERFA
using LinearAlgebra: norm
using SPICE: spkezr, furnsh, kclear
using RemoteFiles: @RemoteFile, @RemoteFileSet, download, path

@RemoteFile(
    DE200,
    "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/a_old_versions/de200.bsp",
    dir=joinpath(@__DIR__(), "data"),
)

download(DE200)
furnsh(path(DE200))

getet(epoch) = value(seconds(j2000(epoch)))

import AstroBase.Ephemerides: position!, velocity!, state!

struct DummyEphemeris <: AbstractEphemeris end
const eph = DummyEphemeris()

function position!(pos, ::DummyEphemeris, ::Epoch, ::NAIFId, ::NAIFId)
    pos .+= [1.0, 2.0, 3.0]
end

function velocity!(vel, ::DummyEphemeris, ::Epoch, ::NAIFId, ::NAIFId)
    vel .+= [1.0, 2.0, 3.0]
end

function state!(pos, vel, ::DummyEphemeris, ::Epoch, ::NAIFId, ::NAIFId)
    (pos .+= [1.0, 2.0, 3.0], vel .+= [1.0, 2.0, 3.0])
end

struct DummyEphemeris2 <: AbstractEphemeris end
const eph2 = DummyEphemeris2()

function position!(pos, ::DummyEphemeris2, ::Epoch, ::NAIFId, ::NAIFId)
    pos .+= [1.0, 2.0, 3.0]
end

function velocity!(vel, ::DummyEphemeris2, ::Epoch, ::NAIFId, ::NAIFId)
    vel .+= [1.0, 2.0, 3.0]
end

@testset "Ephemerides" begin
    @testset "AbstractEphemeris" begin
        ep = TDBEpoch(2000, 1, 1)
        res = [1.0, 2.0, 3.0]
        @test position!(zeros(3), eph, ep, ssb, sun) == res
        @test position(eph, ep, ssb, sun) == res
        @test velocity!(zeros(3), eph, ep, ssb, sun) == res
        @test velocity(eph, ep, ssb, sun) == res
        @test state!(zeros(3), zeros(3), eph, ep, ssb, sun) == (res, res)
        @test state(eph, ep, ssb, sun) == (res, res)

        @test position!(zeros(3), eph, ep, luna, io) == 4res
        @test position(eph, ep, luna, io) == 4res
        @test velocity!(zeros(3), eph, ep, luna, io) == 4res
        @test velocity(eph, ep, luna, io) == 4res
        @test state!(zeros(3), zeros(3), eph, ep, luna, io) == (4res, 4res)
        @test state(eph, ep, luna, io) == (4res, 4res)

        @test position!(zeros(3), eph, ep, sun) == res
        @test position(eph, ep, sun) == res
        @test velocity!(zeros(3), eph, ep, sun) == res
        @test velocity(eph, ep, sun) == res
        @test state!(zeros(3), zeros(3), eph, ep, sun) == (res, res)
        @test state(eph, ep, sun) == (res, res)

        @test position!(zeros(3), eph, ep, ssb) == zeros(3)
        @test position(eph, ep, ssb) == zeros(3)
        @test velocity!(zeros(3), eph, ep, ssb) == zeros(3)
        @test velocity(eph, ep, ssb) == zeros(3)
        @test state!(zeros(3), zeros(3), eph, ep, ssb) == (zeros(3), zeros(3))
        @test state(eph, ep, ssb) == (zeros(3), zeros(3))

        @test position!(zeros(3), eph2, ep, luna, io) == 4res
        @test position(eph2, ep, luna, io) == 4res
        @test velocity!(zeros(3), eph2, ep, luna, io) == 4res
        @test velocity(eph2, ep, luna, io) == 4res
        @test state!(zeros(3), zeros(3), eph2, ep, luna, io) == (4res, 4res)
        @test state(eph2, ep, luna, io) == (4res, 4res)
    end
    @testset "Simon/Bretagnon" begin
        bodies = (
            mercury,
            venus,
            earth_barycenter,
            mars,
            jupiter,
            saturn,
            uranus,
            neptune,
        )
        pos_error = (
            0.5e6,
            1.1e6,
            1.3e6,
            9e6,
            82e6,
            263e6,
            661e6,
            248e6,
        )
        vel_error = (
            0.7,
            0.9,
            1.0,
            2.5,
            8.2,
            24.6,
            27.4,
            21.4,
        )
        ep = TDBEpoch(2009, 5, 24)
        et = getet(ep)
        jd1, jd2 = value.(julian_twopart(ep))
        @testset for (id, body) in enumerate(bodies)
            @testset "ERFA" begin
                r_erfa, v_erfa = ERFA.plan94(jd1, jd2, id)
                r_erfa .*= astronomical_unit()
                v_erfa .*= astronomical_unit()
                v_erfa ./= SECONDS_PER_DAY
                rv_erfa = (r_erfa, v_erfa)
                r_act = position!(zeros(3), simon_bretagnon, ep, sun, body)
                v_act = velocity!(zeros(3), simon_bretagnon, ep, sun, body)
                rv_act = state!(zeros(3), zeros(3), simon_bretagnon, ep, sun, body)
                @testset for i in 1:3
                    @test r_act[i] ≈ r_erfa[i] atol=1.0
                end
                @testset for i in 1:3
                    @test v_act[i] ≈ v_erfa[i] atol=5.0
                end
                @testset for i in 1:3
                    @test rv_act[1][i] ≈ rv_erfa[1][i] atol=1.0
                    @test rv_act[2][i] ≈ rv_erfa[2][i] atol=5.0
                end
            end
            @testset "JPL" begin
                name = string(body)
                if body in (jupiter, saturn, uranus, neptune)
                    name *= " Barycenter"
                end
                s_exp = spkezr(name, et, "J2000", "NONE", "SUN")[1]
                r_exp = s_exp[1:3]
                v_exp = s_exp[4:6]
                rv_exp = (r_exp, v_exp)
                r_act = position!(zeros(3), simon_bretagnon, ep, sun, body)
                v_act = velocity!(zeros(3), simon_bretagnon, ep, sun, body)
                rv_act = state!(zeros(3), zeros(3), simon_bretagnon, ep, sun, body)
                @test norm(r_act) ≈ norm(r_exp) atol=pos_error[id]
                @test norm(v_act) ≈ norm(v_exp) atol=vel_error[id]
                @test norm(rv_act[1]) ≈ norm(rv_exp[1]) atol=pos_error[id]
                @test norm(rv_act[2]) ≈ norm(rv_exp[2]) atol=vel_error[id]
                s_exp = spkezr("MERCURY", et, "J2000", "NONE", "EARTH_BARYCENTER")[1]
                r_exp = s_exp[1:3]
                v_exp = s_exp[4:6]
                rv_exp = (r_exp, v_exp)
                r_act = position!(zeros(3), simon_bretagnon, ep, earth_barycenter, mercury)
                v_act = velocity!(zeros(3), simon_bretagnon, ep, earth_barycenter, mercury)
                rv_act = state!(zeros(3), zeros(3), simon_bretagnon, ep, earth_barycenter, mercury)
                @test norm(r_act) ≈ norm(r_exp) atol=pos_error[3]
                @test norm(v_act) ≈ norm(v_exp) atol=vel_error[3]
                @test norm(rv_act[1]) ≈ norm(rv_exp[1]) atol=pos_error[3]
                @test norm(rv_act[2]) ≈ norm(rv_exp[2]) atol=vel_error[3]
            end
        end
        @test_throws ArgumentError state(simon_bretagnon, ep, earth, luna)
    end
    @testset "VSOP87" begin
        bodies = (
            sun,
            mercury,
            venus,
            earth,
            mars,
            jupiter,
            saturn,
            uranus,
            neptune,
        )
        pos_error = (
            0.5e6,
            0.5e6,
            1.1e6,
            1.3e6,
            9e6,
            82e6,
            263e6,
            661e6,
            248e6,
        )
        vel_error = (
            0.7,
            0.7,
            0.9,
            1.0,
            2.5,
            8.2,
            24.6,
            27.4,
            21.4,
        )
        ep = TDBEpoch(2019, 5, 6)
        et = getet(ep)
        vsop87 = VSOP87()
        @testset for (id, body) in enumerate(bodies)
            name = string(body)
            if body in (jupiter, saturn, uranus, neptune)
                name *= " Barycenter"
            end
            s_exp = spkezr(name, et, "J2000", "NONE", "SSB")[1]
            r_exp = s_exp[1:3]
            v_exp = s_exp[4:6]
            rv_exp = (r_exp, v_exp)
            r_act = position!(zeros(3), vsop87, ep, ssb, body)
            v_act = velocity!(zeros(3), vsop87, ep, ssb, body)
            rv_act = state!(zeros(3), zeros(3), vsop87, ep, ssb, body)
            @testset for i in 1:3
                @test r_act[i] ≈ r_exp[i] atol=pos_error[id]
            end
            @testset for i in 1:3
                @test v_act[i] ≈ v_exp[i] atol=vel_error[id]
            end
            @testset for i in 1:3
                @test rv_act[1][i] ≈ rv_exp[1][i] atol=pos_error[id]
                @test rv_act[2][i] ≈ rv_exp[2][i] atol=vel_error[id]
            end
        end
        s_exp = spkezr("MERCURY", et, "J2000", "NONE", "EARTH")[1]
        r_exp = s_exp[1:3]
        v_exp = s_exp[4:6]
        rv_exp = (r_exp, v_exp)
        r_act = position!(zeros(3), vsop87, ep, earth, mercury)
        v_act = velocity!(zeros(3), vsop87, ep, earth, mercury)
        rv_act = state!(zeros(3), zeros(3), vsop87, ep, earth, mercury)
        @testset for i in 1:3
            @test r_act[i] ≈ r_exp[i] atol=pos_error[3]
        end
        @testset for i in 1:3
            @test v_act[i] ≈ v_exp[i] atol=vel_error[3]
        end
        @testset for i in 1:3
            @test rv_act[1][i] ≈ rv_exp[1][i] atol=pos_error[3]
            @test rv_act[2][i] ≈ rv_exp[2][i] atol=vel_error[3]
        end
        @test_throws ArgumentError state(vsop87, ep, earth, luna)
    end

    kclear()
end


text:
using BinDeps

@BinDeps.setup

# 32-bit
libmetis = library_dependency("libmetis", aliases=["libmetis5"]) # The default library in all the systems is compiled with 32 bit integer and floats

if is_windows()
    using WinRPM
    provides(WinRPM.RPM, "metis", libmetis, os = :Windows)
end

if is_apple()
    if Pkg.installed("Homebrew") === nothing
        error("Homebrew package not installed, please run Pkg.add(\"Homebrew\")")  end
    using Homebrew
    provides(Homebrew.HB, "metis", libmetis, os = :Darwin)
end

provides(AptGet, "libmetis5", libmetis)

provides(Yum, "metis-5.1.0", libmetis)

provides(Sources, URI("http://glaros.dtc.umn.edu/gkhome/fetch/sw/metis/metis-5.1.0.tar.gz"), libmetis)

metisname = "metis-5.1.0"

provides(Sources, URI("http://glaros.dtc.umn.edu/gkhome/fetch/sw/metis/metis-5.1.0.tar.gz"), libmetis)

depsdir = BinDeps.depsdir(libmetis)
prefix = joinpath(depsdir,"usr")
srcdir = joinpath(depsdir,"src",metisname)

provides(SimpleBuild,
         (@build_steps begin
             GetSources(libmetis)
             (@build_steps begin
	         ChangeDirectory(srcdir)
	         (@build_steps begin
	             `make config shared=1 prefix=$prefix`
	             `make`
	             `make install`
                 end)
	     end)
         end), libmetis, os=:Unix)

# 64-bit
# Only done for UNIX
libmetis64 = library_dependency("libmetis64")

provides(Sources, URI("http://glaros.dtc.umn.edu/gkhome/fetch/sw/metis/metis-5.1.0.tar.gz"), libmetis64)
metis64_depsdir = BinDeps.depsdir(libmetis64)
metis64_prefix = joinpath(metis64_depsdir, "usr")
metis64_srcdir = joinpath(metis64_depsdir, "src", metisname)

provides(SimpleBuild,
         (@build_steps begin
             GetSources(libmetis64)
             (@build_steps begin
	         ChangeDirectory(metis64_srcdir)
	         (@build_steps begin
                 `cp ../../metis.h $metis64_srcdir/include/`
	             `make config shared=1 prefix=$metis64_prefix`
	             `make`
	             `make install`
	             `mv ../../usr/lib/libmetis.so ../../usr/lib/libmetis64.so`
                 end)
	     end)
         end), libmetis64, os=:Unix)

@BinDeps.install Dict([(:libmetis => :libmetis), (:libmetis64 => :libmetis64)])