text:
using Test
using GraphletCounting

using LightGraphs

@testset "Degree oriented DAG" begin
    G = SimpleGraph(6)
    add_edge!(G, 1, 2)
    add_edge!(G, 1, 3)
    add_edge!(G, 3, 4)
    add_edge!(G, 4, 5)
    add_edge!(G, 1, 6)
    DG = toDAG(G)

    expectedDG = SimpleDiGraph(6)
    add_edge!(expectedDG, 2, 1)
    add_edge!(expectedDG, 6, 1)
    add_edge!(expectedDG, 3, 1)
    add_edge!(expectedDG, 3, 4)
    add_edge!(expectedDG, 5, 4)

    # # checking adjmats bypasses the sorting of the fadjlist
    @test adjacency_matrix(DG, dir=:out) == adjacency_matrix(expectedDG, dir=:out)
    @test adjacency_matrix(DG, dir=:in) == adjacency_matrix(expectedDG, dir=:in)
end


text:
include("../src/neuro_module.jl")

module ntest
using ..neuro
using PyPlot
PyPlot.pygui(true)
using Random
Random.seed!(10)

function meshgrid(x, y)
    Lx = length(x)
    Ly = length(y)
    T = Float64

    xm = Matrix{T}(undef, Lx, Ly)
    ym = Matrix{T}(undef, Lx, Ly)

    for jj = 1:Ly
        xm[:, jj] = x
    end

    for ii = 1:Lx
        ym[ii, :] = y
    end

    return (xm, ym)
end

function true_fcn(x)
    # y = 0.5 * x[1] + 2.0 * x[2]
    a = 1.0
    b = 0.1
    y = (a - x[1])^2 + b * (x[2] - x[1]^2)^2
    return y
end

function true_fcn_v2(x)
    # y = 0.5 * x[1] + 2.0 * x[2]
    a = 1.0
    b = 0.1
    y1 = (a - x[1])^2 + b * (x[2] - x[1]^2)^2
    y2 = (a - x[1])^2
    return [y1, y2]
end

function create_sample_container()
    x1 = collect( range(-3.0, 4.0; length = 55) )
    x2 = collect( range(-2.0, 2.0; length = 85) )
    (xm, ym) = meshgrid(x1, x2)

    xv = vec(xm)
    yv = vec(ym)

    xx = [xv'; yv']
    yy = vec( mapslices(true_fcn, xx, dims = 1) )
    sample_cont = neuro.MatrixVectorSampleContainer(xx, yy)

    return sample_cont
end

function create_sample_container_v2()
    x1 = collect( range(-3.0, 4.0; length = 55) )
    x2 = collect( range(-2.0, 2.0; length = 55) )
    (xm, ym) = meshgrid(x1, x2)

    xv = vec(xm)
    yv = vec(ym)

    xx = [xv'; yv']
    yy = mapslices(true_fcn_v2, xx, dims = 1)
    sample_cont = neuro.MatrixMatrixSampleContainer(xx, yy)

    return sample_cont
end

#=
n_neurons_per_layers = [2, 30, 30, 1]
act_fcn              = [neuro.TanhActivationFunction(), neuro.TanhActivationFunction(), neuro.IdentityActivationFunction()]
=#

#=
n_neurons_per_layers = [2, 30, 30, 30, 1]
act_fcn              = [
    neuro.TanhActivationFunction(),
    neuro.TanhActivationFunction(),
    neuro.TanhActivationFunction(),
    neuro.IdentityActivationFunction()]
=#

#=
n_neurons_per_layers = [2, 30, 1]
act_fcn              = [
    neuro.TanhActivationFunction(),
    neuro.IdentityActivationFunction()]
=#

n_neurons_per_layers = [2, 10, 2]
act_fcn              = [
    neuro.TanhActivationFunction(),
    neuro.IdentityActivationFunction()]

nn                   = neuro.NeuralNetwork(n_neurons_per_layers, act_fcn)

#=
wm1 = ones(n_neurons_per_layers[2], n_neurons_per_layers[1])
wm2 = ones(n_neurons_per_layers[3], n_neurons_per_layers[2])

b1 = -ones(n_neurons_per_layers[2])
b2 = -ones(n_neurons_per_layers[3])

neuro.copy_in_bias!(nn, b1, 1)
neuro.copy_in_bias!(nn, b2, 2)

neuro.copy_in_weight_matrix!(nn, wm1, 1)
neuro.copy_in_weight_matrix!(nn, wm2, 2)
=#

sample_cont = create_sample_container_v2()

N_sample = length(sample_cont)
p_cv     = 0.8
N        = round(Int, p_cv * N_sample)
neuro.shuffle!(sample_cont)
(sample_cont_cut, cv_cont) = neuro.split(sample_cont, N)

training_options = neuro.TrainingOptions(
    learning_rate = 0.001,
    batch_size    = 30,
    max_iter      = 10_000)

# (nn, iter) = neuro.train!(nn, sample_cont, training_options)

nn = neuro.train2!(nn, sample_cont_cut, training_options)

# perf     = neuro.evaluate_performance(nn, cv_cont)
perf     = neuro.evaluate_performance(nn, sample_cont)
perfstat = neuro.PerformanceStat(perf)

# y_out = neuro.evaluate_container(nn, cv_cont)
y_out = neuro.evaluate_container(nn, sample_cont)

#=
PyPlot.figure()
PyPlot.surf(sample_cont.x[1,:], sample_cont.x[2,:], sample_cont.y, alpha = 0.5)
PyPlot.plot3D(cv_cont.x[1,:], cv_cont.x[2,:], y_out[1,:], marker = :., markersize = 10, linestyle = :none)
=#

PyPlot.figure()
PyPlot.surf(sample_cont.x[1,:], sample_cont.x[2,:], sample_cont.y[1,:], alpha = 0.5)
# PyPlot.plot3D(cv_cont.x[1,:], cv_cont.x[2,:], y_out[1,:], marker = :., markersize = 10, linestyle = :none)
PyPlot.plot3D(sample_cont.x[1,:], sample_cont.x[2,:], y_out[1,:], marker = :., markersize = 10, linestyle = :none)

PyPlot.figure()
PyPlot.surf(sample_cont.x[1,:], sample_cont.x[2,:], sample_cont.y[2,:], alpha = 0.5)
# PyPlot.plot3D(cv_cont.x[1,:], cv_cont.x[2,:], y_out[2,:], marker = :., markersize = 10, linestyle = :none)
PyPlot.plot3D(sample_cont.x[1,:], sample_cont.x[2,:], y_out[2,:], marker = :., markersize = 10, linestyle = :none)

end # module

ext: jl

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
module OpenFOAM
# that is the name of the module, but it's also only different from the
# new type that it defines by case (the type is OpenFoam)

export Lorenz63,OpenFoam,initCase,run,runQ,readMesh,findTimes,readVar,readVarSpec,stringG,rewriteVar,reshapeMesh,writeDict,writeVolScalarField,writeVolVectorField,writeSurfaceScalarField
# println("foamLia: basic openfoam manipulation")

using DataStructures
import Base.run

function stringG(a::Number)
    if float(int(a)) == a
        return string(int(a))
    else
        return string(a)
    end
end

# abstract model type to subclass Lorenz63,OpenFOAM from
abstract Model

type Lorenz63 <: Model
    # x::Array(Float64,1)
    x
    dt::Function
    J::Function
    nVar::Int
    tstep::Float64
    t::Float64
    window::Float64
    # params::Array(Float64,1)
    params
    TLMMethod::String
    # p_f::Array(Float64,1,1)
    pf
    directory::String
end

# functional programming for breakfast
# the outer function returns a function
# which is the lorenz63 model
function outer()
    inner = function(y::Array,t::Float64)
        b,s,r = (8/3,10.0,28.0)
        [s*(y[2]-y[1]),r*y[1]-y[2]-y[1]*y[3],y[1]*y[2]-b*y[3]]
    end
    return inner
end

# this is the lorenz 63 jacobian
J = function(y::Array,t::Float64)
    b,s,r = (8/3,10.0,28.0)
    [-s s 0.0;-y[3]+r -1.0 -y[1];y[2] y[1] -b]
end

# redefine the default constructor
# with some intermediaries
Lorenz63(dt,J,nVar) = Lorenz63([1.0,1.0,1.0],dt,J,nVar,.01,0.0,20.0,[8/3,10.0,28.0],"rk4prime",[],"/Users/andyreagan/work")
# use the above contstruct #dispatching
Lorenz63(J) = Lorenz63(outer(),J,3)
Lorenz63() = Lorenz63(J)

# main (only) type
# want to be able to store all case-related information
type OpenFoam
    caseFolder::String
    controlDict::OrderedDict
    fvSchemes::OrderedDict
    fvSolution::OrderedDict
    setFieldsDict::OrderedDict
    RASProperties::OrderedDict
    transportProperties::OrderedDict
    turbulenceProperties::OrderedDict
    g::OrderedDict
    blockMeshDict::OrderedDict
    T::OrderedDict
    phi::OrderedDict
    U::OrderedDict
    p::OrderedDict
    meshParameters::OrderedDict
    fullMesh::OrderedDict
end

function create_defaultControlDict()
    # default settings
    controlDict = OrderedDict()
    controlDict["application"] = "buoyantBoussinesqPimpleFoam"
    controlDict["startFrom"] = "startTime"
    controlDict["startTime"] = 0
    controlDict["stopAt"] = "endTime"
    controlDict["endTime"] = 100
    controlDict["deltaT"] = .001
    controlDict["writeControl"] = "runTime"
    controlDict["writeInterval"] = 0.1
    controlDict["purgeWrite"] = 0
    controlDict["writeFormat"] = "ascii"
    controlDict["writePrecision"] = 6
    controlDict["writeCompression"] = "off"
    controlDict["timeFormat"] = "general"
    controlDict["timePrecision"] = 6
    controlDict["runTimeModifiable"] = "true"
    controlDict["adjustTime"] = "no"
    controlDict["maxCo"] = 0.5
    controlDict["libs"] = ["\"libOpenFOAM.so\"","\"libsimpleSwakFunctionObjects.so\"","\"libswakFunctionObjects.so\"","\"libgroovyBC.so\""]
    controlDict
end

function create_defaultFvSchemes()
    fvSchemes = OrderedDict()
    fvSchemes["ddtSchemes"] = OrderedDict()
    fvSchemes["ddtSchemes"]["default"] = "Euler"
    # "CrankNicholson 1"

    fvSchemes["gradSchemes"] = OrderedDict()
    fvSchemes["gradSchemes"]["default"] = "Gauss linear"

    fvSchemes["divSchemes"] = OrderedDict()
    fvSchemes["divSchemes"]["default"] = "none"
    fvSchemes["divSchemes"]["div(phi,U)"] = "Gauss upwind"
    fvSchemes["divSchemes"]["div(phi,T)"] = "Gauss upwind"
    fvSchemes["divSchemes"]["div(phi,k)"] = "Gauss upwind"
    fvSchemes["divSchemes"]["div(phi,epsilon)"] = "Gauss upwind"
    fvSchemes["divSchemes"]["div(phi,R)"] = "Gauss upwind"
    fvSchemes["divSchemes"]["div(R)"] = "Gauss linear"
    fvSchemes["divSchemes"]["div((nuEff*dev(T(grad(U)))))"] = "Gauss linear"

    fvSchemes["laplacianSchemes"] = OrderedDict()
    fvSchemes["laplacianSchemes"]["default"] = "none"
    fvSchemes["laplacianSchemes"]["laplacian(nuEff,U)"] = "Gauss linear uncorrected"
    fvSchemes["laplacianSchemes"]["laplacian(Dp,p_rgh)"] = "Gauss linear uncorrected"
    fvSchemes["laplacianSchemes"]["laplacian(alphaEff,T)"] = "Gauss linear uncorrected"
    fvSchemes["laplacianSchemes"]["laplacian(DkEff,k)"] = "Gauss linear uncorrected"
    fvSchemes["laplacianSchemes"]["laplacian(DepsilonEff,epsilon)"] = "Gauss linear uncorrected"
    fvSchemes["laplacianSchemes"]["laplacian(DREff,R)"] = "Gauss linear uncorrected"

    fvSchemes["interpolationSchemes"] = OrderedDict()
    fvSchemes["interpolationSchemes"]["default"] = "linear"

    fvSchemes["snGradSchemes"] = OrderedDict()
    fvSchemes["snGradSchemes"]["default"] = "uncorrected"

    fvSchemes["fluxRequired"] = OrderedDict()
    fvSchemes["fluxRequired"]["default"] = "no"
    fvSchemes["fluxRequired"]["p_rgh"] = ""

    fvSchemes
end

function create_defaultFvSolution()
    fvSolution = OrderedDict()

    fvSolution["solvers"] = OrderedDict()

    fvSolution["solvers"]["p_rgh"] = OrderedDict()
    fvSolution["solvers"]["p_rgh"]["solver"] = "PCG"
    fvSolution["solvers"]["p_rgh"]["preconditioner"] = "DIC"
    fvSolution["solvers"]["p_rgh"]["tolerance"] = "1e-8"
    fvSolution["solvers"]["p_rgh"]["relTol"] = "0.01"

    fvSolution["solvers"]["p_rghFinal"] = OrderedDict()
    fvSolution["solvers"]["p_rghFinal"]["\$p_rgh"] = ""
    fvSolution["solvers"]["p_rghFinal"]["relTol"] = "0"

    fvSolution["solvers"]["\"(U|T|k|epsilon|R)\""] = OrderedDict()
    fvSolution["solvers"]["\"(U|T|k|epsilon|R)\""]["solver"] = "PBiCG"
    fvSolution["solvers"]["\"(U|T|k|epsilon|R)\""]["preconditioner"] = "DILU"
    fvSolution["solvers"]["\"(U|T|k|epsilon|R)\""]["tolerance"] = "1e-6"
    fvSolution["solvers"]["\"(U|T|k|epsilon|R)\""]["relTol"] = "0.1"

    fvSolution["solvers"]["\"(U|T|k|epsilon|R)Final\""] = OrderedDict()
    fvSolution["solvers"]["\"(U|T|k|epsilon|R)Final\""]["\$U"] = ""
    fvSolution["solvers"]["\"(U|T|k|epsilon|R)Final\""]["relTol"] = "0"

    fvSolution["PIMPLE"] = OrderedDict()
    fvSolution["PIMPLE"]["momentumPredicto"] = "no"
    fvSolution["PIMPLE"]["nOuterCorrectors"] = "1"
    fvSolution["PIMPLE"]["nCorrectors"] = "2"
    fvSolution["PIMPLE"]["nNonOrthogonalCorrectors"] = "0"
    fvSolution["PIMPLE"]["pRefCell"] = "0"
    fvSolution["PIMPLE"]["pRefValue"] = "0"

    fvSolution["relaxationFactors"] = OrderedDict()
    fvSolution["relaxationFactors"]["fields"] = OrderedDict()
    fvSolution["relaxationFactors"]["equations"] = OrderedDict()
    fvSolution["relaxationFactors"]["equations"]["\"(U|T|k|epsilon|R)\""] = "1"
    fvSolution["relaxationFactors"]["equations"]["\"(U|T|k|epsilon|R)Final\""] = "1"

    fvSolution
end

function create_defaultSetFieldsDict()
    setFieldsDict = OrderedDict()
    setFieldsDict["defaultFieldValues"] = [""]
    setFieldsDict["regions"] = [""]
    setFieldsDict
end

function create_defaultG()
    g = OrderedDict()
    g["dimensions"] = "[0 1 -2 0 0 0 0]"
    g["value"] = "( 0 0 -9.81)"
    g
end

function create_defaultRASProperties()
    RASProperties = OrderedDict()
    RASProperties["RASModel"] = "\$RASModel"
    RASProperties["turbulence"] = "\$turbulence"
    RASProperties["printCoeffs"] = "on"
    RASProperties
end

function create_defaultTransportProperties()
    transportProperties = OrderedDict()

    # reference numbers:
    # http://www.engineeringtoolbox.com/air-properties-d_156.html
    # http://www.engineeringtoolbox.com/water-thermal-properties-d_162.html

    transportProperties["transportModel"] = "Newtonian"
    transportProperties["Cp0"] = "4.187"
    # Kj/Kg at 300K

    # trying to add some density!
    transportProperties["rho"] = "rho [ 1 -3 0 0 0 0 0 ] 1000"

    # Laminar viscosity
    transportProperties["nu"] = "nu [0 2 -1 0 0 0 0] 1e-06"
    # was 1e-05 for air

    # Thermal expansion coefficient
    transportProperties["beta"] = "beta [0 0 0 -1 0 0 0] 3e-04"
    # 3e-03 for air

    # Reference temperature
    transportProperties["TRef"] = "TRef [0 0 0 1 0 0 0] 300"
    # Kelvin

    # Laminar Prandtl number
    transportProperties["Pr"] = "Pr [0 0 0 0 0 0 0] 5.43"
    # 0.9 for air

    # Turbulent Prandtl number
    transportProperties["Prt"] = "Prt [0 0 0 0 0 0 0] 5.43"
    #assume it's just the same

    transportProperties
end

function create_defaultTurbulenceProperties()
    turbulenceProperties = OrderedDict()
    turbulenceProperties["simulationType"] = "laminar"
    # other options: "RASModel","LESModel"
    turbulenceProperties
end

function create_defaultBlockMeshDict()
    blockMeshDict = OrderedDict()
    blockMeshDict["application"] = "buoyantBoussinesqPimpleFoam"
    blockMeshDict
end

function create_defaultT()
    defaultT = OrderedDict()
    # defaultT["dimensions"] = [0,0,0,1,0,0,0]
    defaultT["dimensions"] = "[0 0 0 1 0 0 0]"

    defaultT["internalField"] = "uniform 300"

    # this should work, but I must not understand deepcopy
    # emptyBC = OrderedDict()
    # emptyBC["type"] = "empty"
    # defaultT["boundaryField"] = OrderedDict()
    # defaultT["boundaryField"]["front"] = deepcopy(emptyBC)
    # defaultT["boundaryField"]["back"] = deepcopy(emptyBC)
    # groovyBC = OrderedDict()
    # groovyBC["type"] = "groovyBC"
    # groovyBC["value"] = "uniform 300"
    # groovyBC["gradientExpression"] = "\"gradT\""
    # groovyBC["fractionExpression"] = "\"0\""
    # groovyBC["variables"] = "\"Text=280;hc=225;gradT=(Text-T)*hc\""
    # groovyBC["timelines"] = ()
    # function setGroovyBC(Text::Int,hc::Int,expression::String)
    #     groovyBC["variables"] = "\"Text=$(Text);hc=$(hc);$(expression)\""
    #     return groovyBC
    # end
    # setGroovyBC(T::Int) = setGroovyBC(T::Int,225,"gradT=(Text-T)*hc")
    # defaultT["boundaryField"]["topinside"] = deepcopy(setGroovyBC(280))
    # println(defaultT)
    # testvar = deepcopy(groovyBC)
    # println(testvar)
    # defaultT["boundaryField"]["topoutside"] = deepcopy(setGroovyBC(280))
    # println(defaultT)
    # println(groovyBC)
    # defaultT["boundaryField"]["bottominside"] = deepcopy(setGroovyBC(340))
    # println(defaultT)
    # println(testvar)
    # println(groovyBC)
    # defaultT["boundaryField"]["bottomoutside"] = deepcopy(setGroovyBC(340))
    # defaultT["boundaryField"]["topoutside"]["variables"] = "\"Text=280;hc=225;gradT=(Text-T)*hc\""
    # println(defaultT)

    defaultT["boundaryField"] = OrderedDict()
    defaultT["boundaryField"]["front"] = OrderedDict()
    defaultT["boundaryField"]["front"]["type"] = "empty"
    defaultT["boundaryField"]["back"] = OrderedDict()
    defaultT["boundaryField"]["back"]["type"] = "empty"
    defaultT["boundaryField"]["topinside"] = OrderedDict()
    defaultT["boundaryField"]["topinside"]["type"] = "groovyBC"
    defaultT["boundaryField"]["topinside"]["value"] = "uniform 300"
    defaultT["boundaryField"]["topinside"]["gradientExpression"] = "\"gradT\""
    defaultT["boundaryField"]["topinside"]["fractionExpression"] = "\"0\""
    defaultT["boundaryField"]["topinside"]["variables"] = "\"Text=280;hc=225;gradT=(Text-T)*hc;\""
    defaultT["boundaryField"]["topinside"]["timelines"] = ()
    defaultT["boundaryField"]["topoutside"] = OrderedDict()
    defaultT["boundaryField"]["topoutside"]["type"] = "groovyBC"
    defaultT["boundaryField"]["topoutside"]["value"] = "uniform 300"
    defaultT["boundaryField"]["topoutside"]["gradientExpression"] = "\"gradT\""
    defaultT["boundaryField"]["topoutside"]["fractionExpression"] = "\"0\""
    defaultT["boundaryField"]["topoutside"]["variables"] = "\"Text=280;hc=225;gradT=(Text-T)*hc;\""
    defaultT["boundaryField"]["topoutside"]["timelines"] = ()
    defaultT["boundaryField"]["bottominside"] = OrderedDict()
    defaultT["boundaryField"]["bottominside"]["type"] = "groovyBC"
    defaultT["boundaryField"]["bottominside"]["value"] = "uniform 300"
    defaultT["boundaryField"]["bottominside"]["gradientExpression"] = "\"gradT\""
    defaultT["boundaryField"]["bottominside"]["fractionExpression"] = "\"0\""
    defaultT["boundaryField"]["bottominside"]["variables"] = "\"Text=340;hc=225;gradT=(Text-T)*hc;\""
    defaultT["boundaryField"]["bottominside"]["timelines"] = ()
    defaultT["boundaryField"]["bottomoutside"] = OrderedDict()
    defaultT["boundaryField"]["bottomoutside"]["type"] = "groovyBC"
    defaultT["boundaryField"]["bottomoutside"]["value"] = "uniform 300"
    defaultT["boundaryField"]["bottomoutside"]["gradientExpression"] = "\"gradT\""
    defaultT["boundaryField"]["bottomoutside"]["fractionExpression"] = "\"0\""
    defaultT["boundaryField"]["bottomoutside"]["variables"] = "\"Text=340;hc=225;gradT=(Text-T)*hc;\""
    defaultT["boundaryField"]["bottomoutside"]["timelines"] = ()

    defaultT
end

function create_defaultPhi()
    defaultPhi = OrderedDict()
    defaultPhi["dimensions"] = "[0 3 -1 0 0 0 0]"

    defaultPhi["internalField"] = "uniform 0"

    defaultPhi["boundaryField"] = OrderedDict()
    defaultPhi["boundaryField"]["front"] = OrderedDict()
    defaultPhi["boundaryField"]["front"]["type"] = "empty"
    defaultPhi["boundaryField"]["front"]["value"] = "uniform 0" # nonuniform 0()
    defaultPhi["boundaryField"]["back"] = OrderedDict()
    defaultPhi["boundaryField"]["back"]["type"] = "empty"
    defaultPhi["boundaryField"]["back"]["value"] = "uniform 0" # nonuniform 0()
    defaultPhi["boundaryField"]["topinside"] = OrderedDict()
    defaultPhi["boundaryField"]["topinside"]["type"] = "calculated"
    defaultPhi["boundaryField"]["topinside"]["value"] = "uniform 0" # nonuniform 0()
    defaultPhi["boundaryField"]["topoutside"] = OrderedDict()
    defaultPhi["boundaryField"]["topoutside"]["type"] = "calculated"
    defaultPhi["boundaryField"]["topoutside"]["value"] = "uniform 0" # nonuniform 0()
    defaultPhi["boundaryField"]["bottominside"] = OrderedDict()
    defaultPhi["boundaryField"]["bottominside"]["type"] = "calculated"
    defaultPhi["boundaryField"]["bottominside"]["value"] = "uniform 0" # nonuniform 0()
    defaultPhi["boundaryField"]["bottomoutside"] = OrderedDict()
    defaultPhi["boundaryField"]["bottomoutside"]["type"] = "calculated"
    defaultPhi["boundaryField"]["bottomoutside"]["value"] = "uniform 0" # nonuniform 0()

    defaultPhi
end

function create_defaultU()
    defaultU = OrderedDict()
    defaultU["dimensions"] = "[0 1 -1 0 0 0 0]"

    defaultU["internalField"] = "uniform (0 0 0)"

    defaultU["boundaryField"] = OrderedDict()
    defaultU["boundaryField"]["front"] = OrderedDict()
    defaultU["boundaryField"]["front"]["type"] = "empty"
    defaultU["boundaryField"]["back"] = OrderedDict()
    defaultU["boundaryField"]["back"]["type"] = "empty"
    defaultU["boundaryField"]["topinside"] = OrderedDict()
    defaultU["boundaryField"]["topinside"]["type"] = "fixedValue"
    defaultU["boundaryField"]["topinside"]["value"] = "uniform (0 0 0)"
    defaultU["boundaryField"]["topoutside"] = OrderedDict()
    defaultU["boundaryField"]["topoutside"]["type"] = "fixedValue"
    defaultU["boundaryField"]["topoutside"]["value"] = "uniform (0 0 0)"
    defaultU["boundaryField"]["bottominside"] = OrderedDict()
    defaultU["boundaryField"]["bottominside"]["type"] = "fixedValue"
    defaultU["boundaryField"]["bottominside"]["value"] = "uniform (0 0 0)"
    defaultU["boundaryField"]["bottomoutside"] = OrderedDict()
    defaultU["boundaryField"]["bottomoutside"]["type"] = "fixedValue"
    defaultU["boundaryField"]["bottomoutside"]["value"] = "uniform (0 0 0)"

    defaultU
end

function create_defaultP()
    defaultP = OrderedDict()
    defaultP["dimensions"] = "[0 2 -2 0 0 0 0]"

    defaultP["internalField"] = "uniform 0"

    defaultP["boundaryField"] = OrderedDict()
    defaultP["boundaryField"]["front"] = OrderedDict()
    defaultP["boundaryField"]["front"]["type"] = "empty"
    defaultP["boundaryField"]["back"] = OrderedDict()
    defaultP["boundaryField"]["back"]["type"] = "empty"
    defaultP["boundaryField"]["topinside"] = OrderedDict()
    defaultP["boundaryField"]["topinside"]["type"] = "calculated"
    defaultP["boundaryField"]["topinside"]["value"] = "uniform 0" # nonuniform 0()
    defaultP["boundaryField"]["topoutside"] = OrderedDict()
    defaultP["boundaryField"]["topoutside"]["type"] = "calculated"
    defaultP["boundaryField"]["topoutside"]["value"] = "uniform 0" # nonuniform 0()
    defaultP["boundaryField"]["bottominside"] = OrderedDict()
    defaultP["boundaryField"]["bottominside"]["type"] = "calculated"
    defaultP["boundaryField"]["bottominside"]["value"] = "uniform 0" # nonuniform 0()
    defaultP["boundaryField"]["bottomoutside"] = OrderedDict()
    defaultP["boundaryField"]["bottomoutside"]["type"] = "calculated"
    defaultP["boundaryField"]["bottomoutside"]["value"] = "uniform 0" # nonuniform 0()

    defaultP
end

function create_defaultMeshParam()
    meshParam = OrderedDict()
    meshParam["baseMeshDir"] = "/users/a/r/areagan/scratch/run/2014-10-23-all-meshes/"
    meshParam["dim"] = 2
    meshParam["x"] = 250
    meshParam["y"] = 40
    meshParam["refinements"] = 0
    meshParam
end

function create_defaultMesh()
    mesh = OrderedDict()
    mesh["points"] = zeros(Float64,3,1) # (-0.0001 -0.375 0)
    mesh["faces"] = zeros(Int64,4,1) # 4(2 84 85 3)
    mesh["owner"] = zeros(Int64,1,1) # 1
    mesh["cellFaces"] = zeros(Int64,6,1) # 1
    mesh["cellCenters"] = zeros(Float64,3,1) # 1
    # mesh["boundary"] = []
    mesh
end

# meshParameters::OrderedDict
# fullMesh::OrderedDict
# this deep copy is not working
# OpenFoam(folder) = OpenFoam(folder,defaultControlDict,defaultTurbulenceProperties,deepcopy(defaultT),defaultMeshParam,defaultMesh)
# OpenFoam(folder) = OpenFoam(folder,defaultControlDict,defaultTurbulenceProperties,create_defaultT(),defaultMeshParam,defaultMesh)
OpenFoam(folder) = OpenFoam(folder,create_defaultControlDict(),create_defaultFvSchemes(),create_defaultFvSolution(),create_defaultSetFieldsDict(),create_defaultRASProperties(),create_defaultTransportProperties(),create_defaultTurbulenceProperties(),create_defaultG(),create_defaultBlockMeshDict(),create_defaultT(),create_defaultPhi(),create_defaultU(),create_defaultP(),create_defaultMeshParam(),create_defaultMesh())

header = """/*--------------------------------*- C++ -*----------------------------------*\
| =========                |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.2.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/"""

lbreak = "// ************************************************************************* //\n\n"

# stolen from dictutils
type CompactRepr
  keySeparator::String
  entrySeparator::String
end

DefaultCompactRepr = CompactRepr(" ", "\n\n")

function showCompact{K,V}(m::OrderedDict{K,V}, config::CompactRepr)
  # Note that without the typehint on the Array, this will somehow get
  # converted to nothing on empty input.
  return join(
    String[string(k) * config.keySeparator * string(v)
           for (k, v) in m],
    config.entrySeparator
  )
end

showCompact{K,V}(m::OrderedDict{K,V}) = showCompact(m, DefaultCompactRepr)
# end stolen from dictutils

function writeDict(o::OpenFoam,d::OrderedDict,name::String,location::String,class::String)
    # mkdir(join([o.caseFolder,location],"/"))

    f = open(join([o.caseFolder,location,name],"/"),"w")

    # write the header
    write(f,string(header,"\n"))

    # write the file info
    finfo = string("FoamFile\n{\n",showCompact(OrderedDict([("version","2.0"),("format","ascii"),("class",class),("location",string("\"",location,"\"")),("object",name)]),CompactRepr(" ", ";\n")),";\n}\n\n")
    write(f,finfo)

    write(f,lbreak)

    # write the main info
    maininfo = string(serializeD(d),"\n\n")
    write(f,maininfo)

    write(f,lbreak)
    close(f)
end

writeDict(o::OpenFoam,d::OrderedDict,name::String,location::String) = writeDict(o::OpenFoam,d::OrderedDict,name::String,location::String,"dictionary")

writeDict(o) = writeDict(o,o.controlDict,"controlDict","system")

function serializeD(d::OrderedDict)
  join(String["\n" * string(k) * " " * serializeDinner(v,1) for (k, v) in d],"\n")
end

function serializeDinner(d::OrderedDict,depth::Int)
  return join(["\n"*repeat("    ",depth-1)*"{\n",join(
    String[repeat("    ",depth) * string(k) * " " * serializeDinner(v,depth+1)
           for (k, v) in d],
    "\n"
  ),"\n"*repeat("    ",depth-1)*"}"],"")
end

function serializeDinner(a::Array,depth::Int)
    if typeof(a[1]) == ASCIIString
        # println("joining array")
        # println(a)
        # println(join(a,"\n"))
        return join(["(",join(a,"\n"),");"],"\n")
    else
        # println("not joining array")
        return join(["[",join(a," "),"];"],"")
    end
end

# function serializeDinner(s::String,depth::Int)
#     join([s,";"],"")
# end

function serializeDinner(s::Any,depth::Int)
    join([string(s),";"],"")
end

function writeVolScalarField(o::OpenFoam,d::OrderedDict,name::String,location::String)
    f = open(join([o.caseFolder,location,name],"/"),"w")

    # write the header
    write(f,string(header,"\n"))

    # write the file info
    finfo = string("FoamFile\n{\n",showCompact(OrderedDict([("version","2.0"),("format","ascii"),("class","volScalarField"),("location",string("\"",location,"\"")),("object",name)]),CompactRepr(" ", ";\n")),";\n}\n\n")
    write(f,finfo)

    write(f,lbreak)

    maininfo = string(serializeD(d),"\n\n")
    write(f,maininfo)

    write(f,lbreak)
    close(f)
end

function writeVolVectorField(o::OpenFoam,d::OrderedDict,name::String,location::String)
    f = open(join([o.caseFolder,location,name],"/"),"w")

    # write the header
    write(f,string(header,"\n"))

    # write the file info
    finfo = string("FoamFile\n{\n",showCompact(OrderedDict([("version","2.0"),("format","ascii"),("class","volVectorField"),("location",string("\"",location,"\"")),("object",name)]),CompactRepr(" ", ";\n")),";\n}\n\n")
    write(f,finfo)

    write(f,lbreak)

    maininfo = string(serializeD(d),"\n\n")
    write(f,maininfo)

    write(f,lbreak)
    close(f)
end

function writeSurfaceScalarField(o::OpenFoam,d::OrderedDict,name::String,location::String)
    f = open(join([o.caseFolder,location,name],"/"),"w")

    # write the header
    write(f,string(header,"\n"))

    # write the file info
    finfo = string("FoamFile\n{\n",showCompact(OrderedDict([("version","2.0"),("format","ascii"),("class","surfaceScalarField"),("location",string("\"",location,"\"")),("object",name)]),CompactRepr(" ", ";\n")),";\n}\n\n")
    write(f,finfo)

    write(f,lbreak)

    maininfo = string(serializeD(d),"\n\n")
    write(f,maininfo)

    write(f,lbreak)
    close(f)
end

# writeDict(o) = writeDict(o,o.controlDict,"controlDict","system")

function copyFromBase(o::OpenFoam,files::Array,baseCase::String)
    for file in files
        s = join([baseCase,file],"/")
        d = join([o.caseFolder,file],"/")
        # println("copying $(s) to $(d)")
        cp(s,d)
	if file == "Allrun"
	    run(`chmod +x $d`)
	end
    end
end

# allFiles = ["Allrun","0/alphat","0/epsilon","0/k","0/nut","0/p","0/p_rgh","0/T","0/U","constant/g","constant/RASProperties","constant/transportProperties","constant/turbulenceProperties","system/controlDict","system/fvSchemes","system/fvSolution","system/FieldsDict"]
allFilesTurbulent = ["Allrun","0/alphat","0/epsilon","0/k","0/nut","0/p","0/p_rgh","0/T","0/U"]
allFiles = ["Allrun","0/alphat","0/p","0/p_rgh","0/T","0/U"]
# meshFiles = [join(["constant","polyMesh",x],"/") for x in ["blockMeshDict","blockMeshDict3D","boundary","faces","neighbour","owner","points"]]
meshFiles = [join(["constant","polyMesh",x],"/") for x in ["boundary","faces","neighbour","owner","points"]]

copyFromBase(o::OpenFoam,baseCase::String) = copyFromBase(o::OpenFoam,append!(allFiles,meshFiles),baseCase::String)

function initCase(o::OpenFoam,baseCase::String)
    # println("making sure folders exist")
    for folder in [join([o.caseFolder,"system"],"/"),join([o.caseFolder,"constant","polyMesh"],"/"),join([o.caseFolder,"0"],"/")]
        if !isdir(folder)
            mkpath(folder)
        else
            println(string(folder," exists"))
        end
    end
    println("copying over base case")
    if o.turbulenceProperties["simulationType"] == "laminar"
        copyFromBase(o,baseCase)
    else
        copyFromBase(o,allFilesTurbulent,baseCase)
        copyFromBase(o,meshFiles,baseCase)
    end

    println("writing out controlDict")
    writeDict(o,o.controlDict,"controlDict","system")

    println("writing out fvSchemes")
    writeDict(o,o.fvSchemes,"fvSchemes","system")

    println("writing out fvSolution")
    writeDict(o,o.fvSolution,"fvSolution","system")

    println("writing out setFieldsDict")
    writeDict(o,o.setFieldsDict,"setFieldsDict","system")

    println("writing out RASProperties")
    writeDict(o,o.RASProperties,"RASProperties","constant")

    println("writing out transportProperties")
    writeDict(o,o.transportProperties,"transportProperties","constant")

    println("writing out turbulenceProperties")
    writeDict(o,o.turbulenceProperties,"turbulenceProperties","constant")

    println("writing out g")
    writeDict(o,o.g,"g","constant","uniformDimensionedVectorField")

    println("writing out T")
    writeVolScalarField(o,o.T,"T","0")

    println("writing out phi")
    writeSurfaceScalarField(o,o.phi,"phi","0")

    println("writing out U")
    writeVolVectorField(o,o.U,"U","0")

    println("writing out p")
    writeVolScalarField(o,o.p,"p","0")
end

# defaultMesh["points"] = [] # (-0.0001 -0.375 0)
# defaultMesh["faces"] = [] # 4(2 84 85 3)
# defaultMesh["neighbour"] = [] # 1
# defaultMesh["owner"] = [] # 0

function readMesh(o::OpenFoam)
    # give this function level scope
    numcells = 0

    # read the mesh from the case
    # goal is to fill up the mesh property

    # first read the points
    f = open(join([o.caseFolder,"constant","polyMesh","points"],"/"),"r")
    b = false
    while !b
        m = match(r"([0-9]+)\n",readline(f))
        if m != nothing
            # println("done with initial read")
            # println("there are $(m.captures[1]) points")
            o.fullMesh["points"] = zeros(3,int(m.captures[1]))
            b = true
        end
    end
    i = 0
    for line in eachline(f)
        m = match(r"\(([-.0-9]+) ([-.0-9]+) ([-.0-9]+)\)\n",line)
        if m != nothing
            i=i+1
            o.fullMesh["points"][:,i] = map(float,m.captures)
        end
    end
    close(f)
    # println("here are some points")
    # println(o.fullMesh["points"][:,100])

    # now read the faces
    f = open(join([o.caseFolder,"constant","polyMesh","faces"],"/"),"r")
    b = false
    while !b
        m = match(r"([0-9]+)\n",readline(f))
        if m != nothing
            # println("done with initial read")
            # println("there are $(m.captures[1]) faces")
            o.fullMesh["faces"] = zeros(Int,4,int(m.captures[1]))
            b = true
        end
    end
    i = 0
    for line in eachline(f)
        m = match(r"4\(([0-9]+) ([0-9]+) ([0-9]+) ([0-9]+)\)\n",line)
        if m != nothing
            i=i+1
            o.fullMesh["faces"][:,i] = map((x)->x+1,map(int,m.captures))
        end
    end
    close(f)
    # println("here is a face's points nums")
    # println(o.fullMesh["faces"][:,100])
    # println("here are those points")
    for p in o.fullMesh["faces"][:,100]
        # println(o.fullMesh["points"][:,p])
    end

    # rather than read into owner, read the faces into their cell
    # each cell needs 6 faces
    f = open(join([o.caseFolder,"constant","polyMesh","owner"],"/"),"r")
    b = false
    while !b
        a = readline(f)
        c = match(r"nCells: ([0-9]+)",a)
        if c != nothing
            # println("there are $(c.captures[1]) cells")
            numcells = int(c.captures[1])
            o.fullMesh["cellFaces"] = zeros(Int,6,numcells)
            o.fullMesh["cellCenters"] = zeros(Float64,3,numcells)
        end
        m = match(r"([0-9]+)\n",a)
        if m != nothing
            # println("done with initial read")
            # println("there are $(m.captures[1]) faces owned")
            o.fullMesh["owner"] = zeros(Int64,int(m.captures[1]))
            b = true
        end
    end
    i = 0
    for line in eachline(f)
        m = match(r"([0-9]+)\n",line)
        if m != nothing
            i=i+1
            j = 1
            while o.fullMesh["cellFaces"][j,int(m.captures[1])+1] != 0
                j=j+1
            end
            o.fullMesh["cellFaces"][j,int(m.captures[1])+1] = i
            o.fullMesh["owner"][i] = int(m.captures[1])+1
        end
    end
    close(f)
    # println(o.fullMesh["cellFaces"][:,100:110])

    f = open(join([o.caseFolder,"constant","polyMesh","neighbour"],"/"),"r")
    b = false
    while !b
        m = match(r"([0-9]+)\n",readline(f))
        if m != nothing
            # println("done with initial read")
            b = true
        end
    end
    i = 0
    for line in eachline(f)
        m = match(r"([0-9]+)\n",line)
        if m != nothing
            i=i+1
            j = 1
            while o.fullMesh["cellFaces"][j,int(m.captures[1])+1] != 0
                j=j+1
            end
            o.fullMesh["cellFaces"][j,int(m.captures[1])+1] = i
        end
    end
    close(f)
    # println(o.fullMesh["cellFaces"][:,100:110])

    for i in 1:numcells
        # add up all of the points, for all of the faces
        # there will be a total of 6*4 = 24 points
        # whereas really we only need 8
        # the triple redundancy won't matter
        # and trying to eliminate it would be slower
        center = zeros(Float64,3)
        for j in 1:6
           for p in o.fullMesh["faces"][:,o.fullMesh["cellFaces"][j,i]]
               center += o.fullMesh["points"][:,p]
           end
        end
        o.fullMesh["cellCenters"][:,i] = center/24
    end

    # println(o.fullMesh["cellCenters"][:,100:110])

    # faces, cells
    # take them while we're at it
    fa,ce = takeSlices(o)
    return fa,ce
end

# here are a couple strings to test of that readVar regex
# teststring = "(2.60429e-31 0.000193236 0.0619409)\n"
# teststring = "(0 8.77659e-05 0.0561036)\n"
# mr = r"\(([0-9.\-e]+) ([0-9.\-e]+) ([0-9.\-e]+)\)\n"
# m = match(mr,teststring)
# if m ≠ nothing
#     print(m)
# else
#     println("you will amount to")
#     print(m)
# end

function readVar(o::OpenFoam,t::String,v::String)
    # just go read that file
    # t is the time (a string!)
    # v is the variable
    cd(o.caseFolder)
    f = open(join([t,v],"/"),"r")
    b = false
    # couple defaults
    n = 1
    mr = r"([0-9]+)\n"
    var = []
    while !b
        a = readline(f)
        c = match(r"class\s+([a-zA-Z]+);",a)
        if c != nothing
            # println("this variable is class $(c.captures[1])")
            if c.captures[1] == "volScalarField"
                mr = r"([0-9.\-e]+)\n"
                n = 1
            elseif c.captures[1] == "surfaceScalarField"
                mr = r"([0-9.\-e]+)\n"
                n = 1
            elseif c.captures[1] == "volVectorField"
                mr = r"\(([0-9.\-e]+) ([0-9.\-e]+) ([0-9.\-e]+)\)\n"
                n = 3
            elseif c.captures[1] == "surfaceVectorField"
                mr = r"\(([0-9.\-e]+) ([0-9.\-e]+) ([0-9.\-e]+)\)\n"
                # mr = r"\(([0-9.\-e]+e*[0-9.\-]+) ([0-9.\-]+e*[0-9.\-]+) ([0-9.\-]+e*[0-9.\-]+)\)\n"
                n = 3
            end
        end
        m = match(r"([0-9]+)\n",a)
        if m != nothing
            # println("done with initial read")
            # println("there are $(m.captures[1]) variables to read")
            var = zeros(Float64,n,int(m.captures[1]))
            b = true
        end
    end
    i = 0
    # println(size(var))
    for line in eachline(f)
        m = match(mr,line)
        if m ≠ nothing && i < size(var)[2]
            i=i+1
            # println(m)
            # println(map(string,m.captures))
            # println(i)
            # println(var[:,i])
            var[:,i] = map(float,map(string,m.captures))
        end
    end
    close(f)
    return var
end

function rewriteVar(o::OpenFoam,t::String,v::String,var::Array)
    # just go read that file
    # t is the time (a string!)
    # v is the variable name, "T"
    # var is the array of the variable

    cd(o.caseFolder)
    f = open(join([t,v],"/"),"r")
    b = false
    header = ""
    # couple defaults
    n = 1
    len = 40000
    mr = r"([0-9]+)\n"
    while !b
        a = readline(f)
        header = join([header,a],"")
        c = match(r"class\s+([a-zA-Z]+);",a)
        if c != nothing
            # println("this variable is class $(c.captures[1])")
            if c.captures[1] == "volScalarField"
                mr = r"([0-9.\-]+e*[0-9.\-]+)\n"
                n = 1
            elseif c.captures[1] == "surfaceScalarField"
                mr = r"([0-9.\-]+e*[0-9.\-]+)\n"
                n = 1
            elseif c.captures[1] == "volVectorField"
                mr = r"\(([0-9.\-]+e*[0-9.\-]+) ([0-9.\-]+e*[0-9.\-]+) ([0-9.\-]+e*[0-9.\-]+)\)\n"
                n = 3
            elseif c.captures[1] == "surfaceVectorField"
                mr = r"\(([0-9.\-]+e*[0-9.\-]+) ([0-9.\-]+e*[0-9.\-]+) ([0-9.\-]+e*[0-9.\-]+)\)\n"
                n = 3
            end
        end
        m = match(r"([0-9]+)\n",a)
        if m != nothing
            # println("done with initial read")
            # println("there are $(m.captures[1]) variables to read")
            len = int(m.captures[1])
            b = true
        end
    end
    i = 0
    println(size(var))

    b = false
    while !b
        a = readline(f)
        if a == ";\n"
            b = true
        end
    end
    # footer = ""
    footer = join(readlines(f),"")
    close(f)
    # println(header)
    # println(footer)

    cp(join([t,v],"/"),join([t,string(v,"_o")],"/"))
    f = open(join([t,v],"/"),"w")
    write(f,header)
    if n == 1
        write(f,"(\n")
        for i in 1:length(var)
            write(f,"$(var[i])\n")
        end
        write(f,")\n")
        write(f,";\n\n")
    elseif n==3
        write(f,"(\n")
        for i in 1:length(var)
            write(f,"($(var[i,1]) $(var[i,2]) $(var[i,3]))\n")
        end
        write(f,")\n")
        write(f,";\n\n")
    end
    write(f,footer)
    close(f)
    return var
end

function readVarSpec(o::OpenFoam,t::String,v::String,ind::Array{Int64})
    # get some specific lines
    # t is the time (a string!)
    # v is the variable
    cd(o.caseFolder)
    f = open(join([t,v],"/"),"r")
    b = false
    # couple defaults
    n = 1
    mr = r"([0-9.\-]+e*[0-9.\-]+)\n"
    var = []
    while !b
        a = readline(f)
        c = match(r"class\s+([a-zA-Z]+);",a)
        if c != nothing
            # println("this variable is class $(c.captures[1])")
            if c.captures[1] == "volScalarField"
                mr = r"([0-9.\-]+e*[0-9.\-]+)\n"
                n = 1
            elseif c.captures[1] == "surfaceScalarField"
                mr = r"([0-9.\-]+e*[0-9.\-]+)\n"
                n = 1
            elseif c.captures[1] == "volVectorField"
                mr = r"\(([0-9.\-]+e*[0-9.\-]+) ([0-9.\-]+e*[0-9.\-]+) ([0-9.\-]+e*[0-9.\-]+)\)\n"
                n = 3
            elseif c.captures[1] == "surfaceVectorField"
                mr = r"\(([0-9.\-]+e*[0-9.\-]+) ([0-9.\-]+e*[0-9.\-]+) ([0-9.\-]+e*[0-9.\-]+)\)\n"
                n = 3
            end
        end
        m = match(r"([0-9]+)\n",a)
        if m != nothing
            # println("done with initial read")
            # println("there are $(m.captures[1]) variables to read")
            var = zeros(Float64,n,length(ind))
            b = true
        end
    end
    i = 0
    j = 1
    while j <= length(ind)
        line = readline(f)
        m = match(mr,line)
        if m != nothing
            i=i+1
            if i == ind[j]
                var[:,j] = map(float,map(string,m.captures))
                j = j+1
            end
        end
    end
    close(f)
    return var
end

function takeSlices(o::OpenFoam)
    # find the points which lie directly on the slices (y or z are 0)
    TOL = 1e-15
    # points = zeros(Int64,4,1) # 12,3,6,9
    # points = [Array(Int64,1) for i in 1:4] # 12,3,6,9
    points = [Int64[] for i in 1:4] # 12,3,6,9
    for i in 1:size(o.fullMesh["points"])[2]
    # for i in 1:length(o.fullMesh["points"][1,:])
        if abs(o.fullMesh["points"][2,i]) < TOL # vertical slice
            if o.fullMesh["points"][3,i] > 0 # top
                append!(points[1],[i])
            else # bottom
                append!(points[3],[i])
            end
        end
        if abs(o.fullMesh["points"][3,i]) < TOL
            if o.fullMesh["points"][2,i] > 0 # right
                append!(points[2],[i])
            else # left
                append!(points[4],[i])
            end
        end
    end
    # return points

    ## make sure that the face contains ALL of these points
    faces = [Int64[] for i in 1:4] # 12,3,6,9
    # go get the faces for which those points are in
    for i in 1:size(o.fullMesh["faces"])[2]
        for j in 1:4
            pointsOnSlice = 0
            for p in o.fullMesh["faces"][:,i]
               if p in points[j]
                   pointsOnSlice = pointsOnSlice+1
                   if pointsOnSlice == 4
                       append!(faces[j],[i])
                   end
                   continue
               else
                   break
               end
            end
        end
    end

    # go get the cells which own those faces (and neighbor them?)
    # really just want the owners...need to go back and save the owner information
    # okay got them
    faceowners = [zeros(Int64,length(faces[1])) for i in 1:4] # 12,3,6,9
    # go get the faces for which those points are in
    for j in 1:4
        for i in 1:length(faces[1])
            faceowners[j][i] = o.fullMesh["owner"][faces[j][i]]
        end
    end
    return faces,faceowners
end

function findTimes(o::OpenFoam)
    times = Int64[]
    cd(o.caseFolder)
    # this works
    # but i don't need to be so redundant
    # for line in split(readall(`ls -R` |> `grep uniform`),"\n")
    #     m = match(r"/([0-9.]+)/uniform",line)
    #     if m != nothing
    #         append!(times,int(m.captures[1]))
    #     end
    # end
    for line in split(readall(`ls`),"\n")
        i = match(r"^([0-9]+)$",line)
        f = match(r"^([0-9]+\.[0-9]+)$",line)
        # println(line)
        # println(m)
        if i != nothing
            # println(i)
            if length(times) > 0
                if typeof(times[1]) == Int64
                    append!(times,[int(i.captures[1])])
                else
                    append!(times,[float(i.captures[1])])
                end
            else
                append!(times,[int(i.captures[1])])
            end
        end
        if f != nothing
            # println(f)
            if length(times) > 0
                if typeof(times[1]) == Int64
                    times = convert(Array{Float64},times)
                end
            else
                times = Float64[]
            end
            append!(times,[float(f.captures[1])])
        end
    end
    sort!(times)
end

qsubheader = """#PBS -l walltime=24:00:00
#PBS -N foamBCTest
#PBS -j oe

cd /users/a/r/areagan/work/2014/11-julia-openfoam

./Allrun"""

function run(o::OpenFoam,c::Cmd,q::String)
    println("submitting the qsub job")
    cd(o.caseFolder)

    walltime = "30:00:00"
    jobname = "foamBCTest"

    f = open("run.qsub","w")

    # write the header
    # write(f,qsubheader)

    a = "#PBS"
    write(f,join([a,"-l",join(["walltime",walltime],"=")]," "))
    write(f,"\n")
    write(f,join([a,"-N",jobname]," "))
    write(f,"\n")
    write(f,join([a,"-j","oe"]," "))
    write(f,"\n\n")
    write(f,"cd $(o.caseFolder)\n\n")
    write(f,"$(string(c)[2:end-1])\n\n")

    close(f)
    run(`qsub -q $(q) run.qsub`)
    # println(`qsub -q $(q) run.qsub`)
end

function run(o::OpenFoam,c::Cmd)
    cd(o.caseFolder)
    run(c)
    0
end

# run function for lorenz model
function run(f::Lorenz63)
    # println("in lorenz function")
    f.x = rk4(f.dt,f.x,[f.t,f.t+f.window],f.tstep)
    f.t += f.window
end

function run(f::Lorenz63,TLM::String)
    # println("in lorenz function, running w TLM")
    f.x,f.pf = rk4p(f.dt,f.J,f.x,[f.t,f.t+f.window],f.tstep)
    f.t += f.window
end

# simple rk4 implementation
function rk4(f::Function,y::Array,tspan::Array,tstep::Float64)
    for t=tspan[1]:tstep:tspan[2]
        s1 = f(y,t)
        s2 = f(y+tstep/2*s1,t+tstep/2)
        s3 = f(y+tstep/2*s2,t+tstep/2)
        s4 = f(y+tstep*s3,t+tstep)
        y += tstep*(s1 + 2*s2 + 2*s3 + s4)/6
    end
    # println(y)
    y
end

# derivative of above method
# also need a function for the jacobian (J here)
function rk4p(f::Function,J::Function,y::Array,tspan::Array,tstep::Float64)
    dim = length(y)
    L = eye(dim)
    for t=tspan[1]:tstep:tspan[2]
        s1 = f(y,t)
        L1 = J(y,t)
        s2 = f(y+tstep/2*s1,t+tstep/2)
        L2 = *(J(y+tstep/2*s1,t+tstep/2),(eye(dim)+tstep/2*L1));
        s3 = f(y+tstep/2*s2,t+tstep/2)
        L3 = *(J(y+tstep/2*s2,t+tstep/2),(eye(dim)+tstep/2*L2));
        s4 = f(y+tstep*s3,t+tstep)
        L4 = *(J(y+tstep*s3,t+tstep),(eye(dim)+tstep*L3));
        L = *(L,eye(dim) + tstep/6*(L1 + 2*L2 + 2*L3 + L4));
        y += tstep*(s1 + 2*s2 + 2*s3 + s4)/6
    end
    y,L
end

function atan3(x,y)
    # like atan2, but go counterclockwise
    # from the left hand side (to always get a positive theta)
    th = atan2(x,y)
    # wrap around
    if th<0
        th = th+2*pi
    end
    th
end

function reshapeMesh(case)
    # don't need the result, but I do want to read the mesh
    readMesh(case)

    x = case.meshParameters["x"]*4
    y = case.meshParameters["y"]+2*case.meshParameters["refinements"]
    points = zeros(Int64,x,y)
    indices = zeros(Int64,x*y,2)

    # println("the mesh is $(x) by $(y)")

    theta = [y for y in 2*pi/x/2:2*pi/x:2*pi-2*pi/x/2]

    # println(size(theta))
    # println(theta[1:10])
    # println(theta[end-10:end])

    # println(size(case.fullMesh["cellCenters"]))
    # println(case.fullMesh["cellCenters"][:,1])

    TOL = 1e-3
    # println(size(case.fullMesh["cellCenters"])[2])
    # println(case.fullMesh["cellCenters"][:,size(case.fullMesh["cellCenters"])[2]])
    for i in 1:size(case.fullMesh["cellCenters"])[2]
        # println("i is $(i)")
        # start them at the right side, go counter clockwise
        # arctan(z/y) where y is right the right (adjacent),
        # z is up (opposite)
        th = atan3(case.fullMesh["cellCenters"][3,i],case.fullMesh["cellCenters"][2,i])

        r = sqrt(case.fullMesh["cellCenters"][3,i]^2+case.fullMesh["cellCenters"][2,i]^2)
        # println("th is $(th)")
        # println("r is $(r)")
        # find the index of theta
        j = 1
        # while abs(th-theta[j])>abs(th-theta[j+1]) > TOL
        while abs(th-theta[j]) > abs(th-theta[j+1]) # && j <= 999
            # println(j)
            # println(abs(th-theta[j]))
            j += 1
            # this can't be the best way
            if j == 1000
                break
            end
        end
        # println("j is $(j)")
        # println(points[j,:])
        k = 1
        while points[j,k] != 0.0
            k += 1
        end
        # println("k is $(k)")
        points[j,k] = i
        # (no point in filling indices until sorted)
    end

    # now go and sort the points by their r value
    for i in 1:x
        ps = points[i,:]
        # go get the r for each of the points
        rs = [sqrt(case.fullMesh["cellCenters"][3,p]^2+case.fullMesh["cellCenters"][2,p]^2) for p in ps]
        perm = sortperm(rs)
        points[i,:] = ps[perm]
        j = 1
        for p in ps
            indices[p,:] = [i,j]
            j+=1
        end
    end

    points,indices
end

end


text:
# Two-qubit stabilizer states from Table II in
# Garcia et al. Efficient Inner-product Algorithm for Stabilizer States.
# arXiv:1210.6646 (2012)
TWOQUBIT = Ket.(normalize!.(Vector{Complex{Float64}}[
    # Separable states
    [0,0,1,0],[0,0,0,1],[0,0,1,1],[0,0,1,-1],[0,0,1,1im],[0,0,1,-1im],
    [1,0,-1,0],[0,1,0,-1],[1,1,-1,-1],[1,-1,-1,1],[1,1im,-1,-1im],[1,-1im,-1,1im],
    [1,0,-1im,0],[0,1,0,-1im],[1,1,-1im,-1im],[1,-1,-1im,1im],[1,1im,-1im,1],[1,-1im,-1im,-1],
    [1,0,1,0],[0,1,0,1],[1,1,1,1],[1,-1,1,-1],[1,1im,1,1im],[1,-1im,1,-1im],
    [1,0,1im,0],[0,1,0,1im],[1,1,1im,1im],[1,-1,1im,-1im],[1,1im,1im,-1],[1,-1im,1im,1],
    [1,0,0,0],[0,1,0,0],[1,1,0,0],[1,-1,0,0],[1,1im,0,0],[1,-1im,0,0],
    # Entangled states
    [0,1,1,0],[1,0,0,-1],[1,0,0,1],[0,1,-1,0],
    [1,0,0,1im],[0,1,1im,0],[0,1,-1im,0],[1,0,0,-1im],
    [1,1,1,-1],[1,1,-1,1],[1,-1,1,1],[1,-1,-1,-1],
    [1,1im,1,-1im],[1,1im,-1,1im],[1,-1im,1,1im],[1,-1im,-1,-1im],
    [1,1,1im,-1im],[1,1,-1im,1im],[1,-1,1im,1im],[1,-1,-1im,-1im],
    [1,1im,1im,1],[1,1im,-1im,-1],[1,-1im,1im,-1],[1,-1im,-1im,1],
    ]),((2,2),))


text:
using Convex
using SCS

# for submodular functions of the form
# f = g(|X|)

a = Variable()
b = Variable()
c = Variable()

g = -3*a + 3*b - c

monotonicity = [a >= 0; b >= a; c >= b]
submodularity = [b - a <= a; c - b <= b - a]
scaling = [a == 1]

println("Maximized redundancy")
problem = maximize(g, [monotonicity; submodularity; scaling])
solver = SCSSolver(verbose = 0)
solve!(problem, solver)

println("Opt=$(problem.optval)")
println("a=$(evaluate(a))")
println("b=$(evaluate(b))")
println("c=$(evaluate(c))")

println()
println("Minimized redundancy")
problem = minimize(g, [monotonicity; submodularity; scaling])
solver = SCSSolver(verbose = 0)
solve!(problem, solver)
println("Opt=$(problem.optval)")
println("a=$(evaluate(a))")
println("b=$(evaluate(b))")
println("c=$(evaluate(c))")

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
import DFControl: cell, atoms, data

function calc_observables(structure::WanStructure{T}, kpoints::Vector{Vec3{T}}, soc::Bool=false) where T
    nat = length(atoms(structure))
    matdim = getwandim(structure)
    klen = length(kpoints)
    calc_angmoms!(structure)
    Sx, Sy, Sz = calc_spins(structure)
    outbands = soc ? wannierbands(2*matdim, kpoints) : wannierbands(matdim, kpoints)
    Threads.@threads for i=1:klen
        k = kpoints[i]
        t_hami, dips = hami_dip_from_k(structure.tbhamis[1], structure.tbRmns[1], k)

        hami = soc ? construct_soc_hami(t_hami, structure) : t_hami
        eigvals, eigvecs = sorted_eig(hami)
        eigvals_k = real(eigvals)
        cm_k      = eigcm(dips, eigvecs)
        L_k, S_k  = eigangmomspin(eigvecs, atoms(structure), Sx, Sy, Sz)
        for j = 1:length(eigvals)
            ob = outbands[j]
            ob.eigvals[i] = eigvals_k[j]
            ob.eigvec[i]  = eigvecs[:, j]
            ob.angmoms[i] = [l[j] for l in L_k]
            ob.cms[i]     = cm_k[j]
            ob.spins[i]   = [s[j] for s in S_k]
        end
    end

    return outbands
end

function calc_observables(structure::WanStructure{T}, kpoints, k_range::StepRangeLen, args...) where T
    mid = div(size(kpoints)[1],2)+1
    beg = Int64(k_range[1])
    steps= div(size(k_range)[1],2)
    last = Int64(k_range[end])
    kxs = [range(kpoints[beg][1], stop = kpoints[mid][1], length = steps)..., kpoints[mid][1], range(kpoints[mid][1], stop = kpoints[last][1], length = steps)[2:end]...]
    kys = [range(kpoints[beg][2],stop = kpoints[mid][2], length = steps)...,kpoints[mid][2] ,range(kpoints[mid][2], stop = kpoints[last][2], length = steps)[2:end]...]
    kzs = [range(kpoints[beg][3],stop = kpoints[mid][3], length = steps)...,kpoints[mid][3] ,range(kpoints[mid][3],stop = kpoints[last][3],length = steps)[2:end]...]
    kxs_t = [range(kpoints[beg][1],stop = kpoints[mid][1],length = steps*100)... ,range(kpoints[mid][1],stop = kpoints[last][1],length = steps*100)[2:end]...]
    kys_t = [range(kpoints[beg][2],stop = kpoints[mid][2],length = steps*100)... ,range(kpoints[mid][2],stop = kpoints[last][2],length = steps*100)[2:end]...]
    kzs_t = [range(kpoints[beg][3],stop = kpoints[mid][3],length = steps*100)... ,range(kpoints[mid][3],stop = kpoints[last][3],length = steps*100)[2:end]...]
    # kxs[div(length(kxs),2)]+=0.00001

    kxs[div(length(kxs),2)] = kxs_t[div(length(kxs_t),2)]
    kxs[div(length(kxs),2)+1] = kxs_t[div(length(kxs_t),2)+2]
    kys[div(length(kxs),2)] = kys_t[div(length(kxs_t),2)]
    kys[div(length(kxs),2)+1] = kys_t[div(length(kxs_t),2)+2]
    kzs[div(length(kxs),2)] = kzs_t[div(length(kxs_t),2)]
    kzs[div(length(kxs),2)+1] = kzs_t[div(length(kxs_t),2)+2]
    kpoints = [[kxs[i],kys[i],kzs[i]] for i=1:length(kxs)]
    return calc_observables(structure, Vec3{T}.(kpoints), args...)
end

function hami_dip_from_k(tbhami, tbdip, k::Vec3{T}) where T
    matdim = size(tbhami[1].block)
    outham = zeros(Complex{T}, matdim)
    outdip = zeros(Point3{Complex{T}}, matdim)
    for i = 1:length(tbhami)
        R_cryst = tbhami[i].R_cryst
        hb = tbhami[i].block
        db = tbdip[i].block
        factor =  ℯ^(-2im*pi*(R_cryst ⋅ k))
        outham .+= factor .* hb
        outdip .+= factor .* db
    end
    outdip_ = real(outdip)
    return outham, [outdip_ fill(zero(Point3{T}), matdim);fill(zero(Point3{T}), matdim) outdip_]
end

function eigangmomspin(eigvecs, atoms::Vector{Atom}, Sx, Sy, Sz)
    outL = Vector{Vec3{Float64}}[]
    outS = Vector{Vec3{Float64}}[]
    len = size(eigvecs)[1]
    for (a, at) in enumerate(atoms)
        L_t    = Vec3{Float64}[]
        S_t    = Vec3{Float64}[]
        for i = 1:len
            eigvec = eigvecs[:, i]
            L = zero(Vec3{Complex{Float64}})
            S = zero(Vec3{Complex{Float64}})
            for i1 = 1:len
                at_index1, blockid1, wfi1 = wfc2atindex(atoms, i1)
                at_index1 != a && continue
                c1 = eigvec[i1]
                for i2 = 1:len
                    at_index2, blockid2, wfi2= wfc2atindex(atoms, i2)
                    at_index2 != a && continue
                    at = atoms[at_index1]
                    c2 = eigvec[i2]

                    if blockid1 == blockid2
                        L += conj(c1) * c2 * angmom(at)[wfi1, wfi2]::Vec3{Complex{Float64}}
                    else
                        S += conj(c1) * c2 * Vec3(Sx[i1, i2], Sy[i1, i2], Sz[i1, i2])
                    end
                end
            end
            push!(L_t, real(L))
            push!(S_t, real(S))
        end
        push!(outL, L_t)
        push!(outS, S_t)
    end
    return outL, outS
end

function eigcm(dip::AbstractMatrix{Point3{T}}, eigvecs) where T
    out = Point3{T}[]
    for i = 1:size(eigvecs)[1]
        eigvec = eigvecs[:, i]
        en = enumerate(eigvec)
        d = zero(Point3{Complex{T}})
        for (i1, v1) in en, (i2, v2) in en
            d += conj(v1) * v2 * dip[i1, i2]
        end
        push!(out, real(d))
    end
    return out
end

function wfc2atindex(atoms::Vector{Atom}, wfci)
    wfcounter = 0
    i = 1
    outi = 1
    while true
        len = length(wfcs(atoms[i]))
        wfcounter += len
        if wfcounter >= wfci
            return i, div1(outi, length(atoms)), len - wfcounter + wfci
        else
            outi += 1
            i +=1
            i = mod1(i, length(atoms))
        end
    end
end


text:
# this is the old upper triangular code

function one_step(A::SMat{T}, sr = 1) where T
  i = sr
  assert(i>0)
  all_r = Array{Int}(0)
  min = A.c
  while i <= length(A.rows)
    @assert length(A.rows[i].pos)>0
    @assert A.rows[i].pos[1] >= sr
    if A.rows[i].pos[1] < min
      min = A.rows[i].pos[1]
      all_r = [i]
    elseif A.rows[i].pos[1] == min
      push!(all_r, i)
    end
    i += 1
  end

  if length(all_r) == 0
    return A.r+1
  end
#  println("found multiple $(length(all_r)) possibilities of weight $([length(A.rows[_i].pos) for _i in all_r])")
  #sorting
  sort!(all_r, lt=function(a,b) return length(A.rows[a].pos) < length(A.rows[b].pos)end)
#  println("found multiple $(length(all_r)) possibilities of weight $([length(A.rows[_i].pos) for _i in all_r])")
  # we need to have rows[all_r[1]] == sr, so that the new pivot will be in
  # row sr

  if !(sr in all_r)
    q = A.rows[sr]
    A.rows[sr] = A.rows[all_r[1]]
    A.rows[all_r[1]] = q
    all_r[1] = sr
  else
    p = findfirst(all_r, sr)
    if p > 1
      q = A.rows[sr]
      A.rows[sr] = A.rows[all_r[1]]
      A.rows[all_r[1]] = q
      all_r[p] = all_r[1]
      all_r[1] = sr
    end
  end
  if length(all_r) == 1
    return sr+1
  end

  for j=2:length(all_r)
    x = A.rows[sr].values[1]
    y = A.rows[all_r[j]].values[1]
    g = x
    @assert x!= 0
    if y % x == 0
      c = -div(y, x)
#      @assert x*c == -y
      add_scaled_row!(A, sr, all_r[j], c)
      @assert !iszero(A[sr])
    else
      g,a,b = gcdx(x, y)
#      @assert g == a*x + b*y
      c = -div(y, g)
#      @assert y == -c*g
      d = div(x, g)
#      @assert x == d*g
      transform_row!(A, sr, all_r[j], a, b, c, d)
      @assert !iszero(A[sr])
    end

#    @assert A.rows[sr].entry[1]valuesval == g
#    @assert A.rows[sr].entry[1].val != 0
#    if length(A.rows[all_r[j]].entry) == 0 ||
#             A.rows[all_r[j]].entry[1].col > min
#    else
#      println("bad", x, " and ", y, " in ", sr, ", ", all_r[j])
#      println("min: ", min)
#    end

#    @assert length(A.rows[all_r[j]].entry) == 0 ||
#             A.rows[all_r[j]].entry[1].col > min
#  println("in one step: ilog2(max) now ", nbits(maximum(abs, A)), " j:", j, " length: ", length(all_r))
  end
  sort!(all_r)
  for j=length(all_r):-1:2
    if iszero(A.rows[all_r[j]])
      deleteat!(A.rows, all_r[j])
      A.r -= 1
    end
  end
  return sr+1
end

doc"""
  upper_triangular{T}(A::SMat{T}; mod = 0)

  Inplace: transform A into an upper triangular matrix. If mod
  is non-zero, reduce entries modulo mod during the computation.
"""
function upper_triangular(A::SMat{T}; mod = 0) where T
  for i = 1:min(rows(A), cols(A))
    x = one_step(A, i)
#    println("after one step: ilog2(max) now ", nbits(maximum(abs, A)))
    if x>A.r
      return
    end
    if A.nnz > (A.r-i) * (A.c-i) /2 || nbits(maximum(abs, A)) > 200
      #println("calling  at level ", i, " bits: ", nbits(maximum(abs, A)), "nnz: ", A.nnz)
      h = sub(A, i:A.r, i:A.c)
      deleteat!(A.rows, i:A.r)
      A.r -= length(i:A.r)
      @assert length(A.rows) == A.r
      h = fmpz_mat(h)
      #println("calling dense hnf on a ", rows(h), " by ", cols(h), " matrix")
      if mod==0
        h = hnf(h)
      else
        h = nmod_mat(mod, h)
        rref!(h)
        h = Array(lift(h))
      end
      h = SMat(h)
      for j in h.rows
        rw = SRow{T}()
        for e in 1:length(j.pos)
          push!(rw.pos, j.pos[e] + i-1)
          push!(rw.values, j.values[e])
        end
        push!(A.rows, rw)
        A.r += 1
      end

      return
    end
  end
  return
end

################################################################################
#
#  One echelonization step w/o transformations
#
################################################################################

# Perform one step in the echelonization, that is, find a pivot and clear
# elements below.
# sr = starting row
function _one_step(A::SMat{T}, sr = 1) where T
  nr, trafo = _one_step_with_trafo(A, sr)
  return nr
end

function _one_step_with_trafo(A::SMat{T}, sr = 1) where T
  trafos = Trafo[]
  i = sr
  assert(i>0)
  all_r = Array{Int}(0)
  min = A.c

  while i <= length(A.rows)
    @assert length(A.rows[i].pos)>0
    @assert A.rows[i].pos[1] >= sr
    if A.rows[i].pos[1] < min
      min = A.rows[i].pos[1]
      all_r = [i]
    elseif A.rows[i].pos[1] == min
      push!(all_r, i)
    end
    i += 1
  end

  if length(all_r) == 0
    return A.r+1, trafos
  end

  #sorting
  sort!(all_r, lt = (a,b) -> (length(A.rows[a].pos) < length(A.rows[b].pos)))
  # we need to have rows[all_r[1]] == sr, so that the new pivot will be in
  # row sr

  if !(sr in all_r)
    push!(trafos, TrafoSwap{T}(sr, all_r[1]))
    q = A.rows[sr]
    A.rows[sr] = A.rows[all_r[1]]
    A.rows[all_r[1]] = q
    all_r[1] = sr
  else
    p = findfirst(all_r, sr)
    if p > 1
      push!(trafos, TrafoSwap{T}(sr, all_r[1]))

      q = A.rows[sr]
      A.rows[sr] = A.rows[all_r[1]]
      A.rows[all_r[1]] = q

      all_r[p] = all_r[1]
      all_r[1] = sr
    end
  end

  if length(all_r) == 1
    return sr + 1, trafos
  end

  for j = 2:length(all_r)
    x = A.rows[sr].values[1]
    y = A.rows[all_r[j]].values[1]
    t = echelonize_basecase!(x, y, sr, all_r[j], A)
    push!(trafos, t)
  end

  sort!(all_r)

  for j = length(all_r):-1:2
    if length(A.rows[all_r[j]].pos) == 0
      deleteat!(A.rows, all_r[j])
      push!(trafos, TrafoDeleteZero{T}(all_r[j]))
      A.r -= 1
    end
  end
  return sr + 1, trafos
end

################################################################################
#
#  Echelonization via dense matrix
#
################################################################################

# The echelonization_via_dense must not remove zero rows automatically
function echelonize_via_dense(h::SMat{nmod})
  R = parent(h.rows[1].values[1])
  # turn into dense structure
  hdense = nmod_mat(h)
  # echelonize
  # I need to the following trick since we do not have the transformation
  rref!(hdense)
  # put back into sparse matrix over R (!= base_ring h)
  return _SMat(hdense, R = R)
end

function echelonize_via_dense_with_trafo(h::SMat{fmpz})
  hdense = fmpz_mat(h)
  # echelonize
  hdense, U = hnf_with_transform(hdense)
  # put back into sparse matrix
  h = _SMat(hdense, zerorows = true)
  return h, U
end

function echelonize_via_dense(h::SMat{fmpz})
  hdense = fmpz_mat(h)
  # echelonize
  hdense = hnf(hdense)
  # put back into sparse matrix
  h = _SMat(hdense)
  return h
end

################################################################################
#
#  Base case of echelonization
#
################################################################################

# Given two rows i, j with entries x, y (x != 0) in the same column, this
# function must produce a zero at entry y.

function echelonize_basecase!(x::T, y::T, i::Int, j::Int, A::Hecke.SMat{T}) where T
  c = -divexact(y, x)
  add_scaled_row!(A, i, j, c)
  return TrafoAddScaled(i, j, c)
end

function echelonize_basecase!(x::fmpz, y::fmpz, i::Int, j::Int,
                              A::Hecke.SMat{fmpz})
  if y % x == 0
    c = -div(y, x)
    add_scaled_row!(A, i, j, c)
    return TrafoAddScaled(i, j, c)
  else
    g,a,b = gcdx(x, y)
    c = -div(y, g)
    d = div(x, g)
    Hecke.transform_row!(A, i, j, a, b, c, d)
    return TrafoParaAddScaled(i, j, a, b, c, d)
  end
end

################################################################################
#
#  Echelonization for sparse matrices
#
################################################################################

# there is some density limit at level i, at which dense echelonization
# will be invoked
function upper_triangular_with_trafo!(M::SMat, density_limit::Float64 = 0.5)
  f = (A, i) -> (A.nnz > (A.r-i) * (A.c-i) * density_limit)
  return _upper_triangular_with_trafo!(M, f)
end

function upper_triangular_with_trafo!(M::SMat{fmpz},
                                      density_limit::Float64 = 0.5,
                                      size_limit::Int = 200)
  f =  (A, i) -> (A.nnz > (A.r-i) * (A.c-i) * density_limit || nbits(maximum(abs, A)) > size_limit)
  return _upper_triangular_with_trafo!(M, f)
end

# is_dense_enough is a function (::SMat{T}, i::Int) -> Bool
# At each level i, is_dense_enough(A, i) is called.
# If it evaluates to true, then dense echelonization will be called.
function _upper_triangular_with_trafo!(A::SMat{T}, is_dense_enough::Function) where T
  trafo = Trafo[]

  for i = 1:min(rows(A), cols(A))
    x, t = _one_step_with_trafo(A, i)
    append!(trafo, t)

    if x > A.r
      return trafo
    end

    if is_dense_enough(A, i)
      h = sub(A, i:A.r, i:A.c)

      h, U = echelonize_via_dense_with_trafo(h)

      # h will have zero rows

      push!(trafo, TrafoPartialDense(i, i:A.r, i:A.c, U))

      deleteat!(A.rows, i:A.r)
      A.r -= length(i:A.r)
      @assert length(A.rows) == A.r

      # append all rows of h to A
      # including zero rows

      for k in 1:length(h.rows)
        j = h.rows[k]
        if length(j.pos) == 0
          push!(A.rows, SRow{T}())
          A.r += 1
          continue
        end
        rw = SRow{T}()
        for e in 1:length(j.pos)
          push!(rw.pos, j.pos[e] + i - 1)
          push!(rw.values, j.values[e])
        end
        push!(A.rows, rw)
        A.r += 1
      end

      k = A.r

      # Now delete trailing zero rows.
      while length(A.rows[k].pos) == 0
        deleteat!(A.rows, k)
        A.r -= 1
        push!(trafo, TrafoDeleteZero{T}(k))
        k -= 1
      end

      # I guess I can do something more clever then the following:
      A.nnz = 0
      for r in A.rows
        A.nnz += length(r.pos)
      end

      return trafo
    end
  end
  return trafo
end

# now the same without transformation

function upper_triangular!(M::SMat, density_limit::Float64 = 0.5)
  f = (A, i) -> (A.nnz > (A.r-i) * (A.c-i) * density_limit)
  return _upper_triangular!(M, f)
end

function upper_triangular!(M::SMat{fmpz}, density_limit::Float64 = 0.5, size_limit::Int = 200)
  f =  (A, i) -> (A.nnz > (A.r-i) * (A.c-i) * density_limit || nbits(maximum(abs, A)) > size_limit)
  return _upper_triangular!(M, f)
end

function _upper_triangular!(A::SMat{T}, is_dense_enough) where T
  for i = 1:min(rows(A), cols(A))
    x = _one_step(A, i)

    if x>A.r
      return nothing
    end

    if is_dense_enough(A, i)

      h = sub(A, i:A.r, i:A.c)
      deleteat!(A.rows, i:A.r)
      A.r -= length(i:A.r)
      @assert length(A.rows) == A.r

      # make dense matrix
      h = echelonize_via_dense(h)

      for j in h.rows
        rw = SRow{T}()
        for e in 1:length(j.pos)
          push!(rw.pos, j.pos[e] + i - 1)
          push!(rw.values, j.values[e])
        end
        push!(A.rows, rw)
        A.r += 1
      end

      # I guess I can do something more clever then the following:
      A.nnz = 0
      for r in A.rows
        A.nnz += length(r.pos)
      end

      return nothing
    end
  end
  return nothing
end

################################################################################
#
#   Smith normal form
#
################################################################################

function _snf_upper_triangular_with_trafo(A::SMat{fmpz})
  # We assume that A is quadratic and upper triangular
  essential_index = 1

  for i in 1:rows(A)
    @assert A.rows[i].pos[1] == i
    if A.rows[i].values[1] != 1
      essential_index = i
      break
    end
  end

  # Let us pick up the transformations

  trafos_right = []
  trafos_left = []

  for i in 1:essential_index-1
    for (j, a) in zip(A.rows[i].pos, A.rows[i].values)
      if j == i
        continue
      end
      push!(trafos_right, TrafoAddScaled(j, i, -a))
    end
  end

  essential_part = fmpz_mat(sub(A, essential_index:rows(A), essential_index:cols(A)))

  snfofess, ltr, rtr = snf_with_transform(essential_part, true, true)

  push!(trafos_left, TrafoPartialDense(essential_index, essential_index:rows(A), essential_index:cols(A), ltr))
  push!(trafos_right, TrafoPartialDense(essential_index, essential_index:rows(A), essential_index:cols(A), rtr))

  return snfofess, trafos_left, trafos_right
end

doc"""
    elementary_divisors(A::SMat{fmpz}) -> Array{fmpz, 1}

> The elementary divisors of $A$, ie. the diagonal elements of the Smith normal
> form of $A$. $A$ needs to be of full rank - currently.
"""
function elementary_divisors(A::SMat{fmpz})
  A = copy(A)
  upper_triangular(A)

  essential_index = 1

  for i in 1:rows(A)
    @assert A.rows[i].pos[1] == i
    if A.rows[i].values[1] != 1
      essential_index = i
      break
    end
  end

  essential_part = fmpz_mat(sub(A, essential_index:rows(A), essential_index:cols(A)))

  s = snf(essential_part)

  return vcat(fmpz[1 for i=1:essential_index-1], fmpz[s[i,i] for i=1:rows(s)])
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

"""
    SystemError(prefix::AbstractString, [errno::Int32])

A system call failed with an error code (in the `errno` global variable).
"""
type SystemError <: Exception
    prefix::AbstractString
    errnum::Int32
    extrainfo
    SystemError(p::AbstractString, e::Integer, extrainfo) = new(p, e, extrainfo)
    SystemError(p::AbstractString, e::Integer) = new(p, e, nothing)
    SystemError(p::AbstractString) = new(p, Libc.errno())
end

"""
    ParseError(msg)

The expression passed to the `parse` function could not be interpreted as a valid Julia
expression.
"""
type ParseError <: Exception
    msg::AbstractString
end

"""
    ArgumentError(msg)

The parameters to a function call do not match a valid signature. Argument `msg` is a
descriptive error string.
"""
type ArgumentError <: Exception
    msg::AbstractString
end

#type UnboundError <: Exception
#    var::Symbol
#end

"""
    KeyError(key)

An indexing operation into an `Associative` (`Dict`) or `Set` like object tried to access or
delete a non-existent element.
"""
type KeyError <: Exception
    key
end

"""
    MethodError(f, args)

A method with the required type signature does not exist in the given generic function.
Alternatively, there is no unique most-specific method.
"""
type MethodError <: Exception
    f
    args
    world::UInt
    MethodError(f::ANY, args::ANY, world::UInt) = new(f, args, world)
end
MethodError(f::ANY, args::ANY) = MethodError(f, args, typemax(UInt))

"""
    EOFError()

No more data was available to read from a file or stream.
"""
type EOFError <: Exception end

"""
    DimensionMismatch([msg])

The objects called do not have matching dimensionality. Optional argument `msg` is a
descriptive error string.
"""
type DimensionMismatch <: Exception
    msg::AbstractString
end
DimensionMismatch() = DimensionMismatch("")

"""
    AssertionError([msg])

The asserted condition did not evaluate to `true`.
Optional argument `msg` is a descriptive error string.
"""
type AssertionError <: Exception
    msg::AbstractString
    AssertionError() = new("")
    AssertionError(msg) = new(msg)
end

#Generic wrapping of arbitrary exceptions
#Subtypes should put the exception in an 'error' field
abstract WrappedException <: Exception

"""
    LoadError(file::AbstractString, line::Int, error)

An error occurred while `include`ing, `require`ing, or `using` a file. The error specifics
should be available in the `.error` field.
"""
type LoadError <: WrappedException
    file::AbstractString
    line::Int
    error
end

"""
    InitError(mod::Symbol, error)

An error occurred when running a module's `__init__` function. The actual error thrown is
available in the `.error` field.
"""
type InitError <: WrappedException
    mod::Symbol
    error
end

ccall(:jl_get_system_hooks, Void, ())

==(w::WeakRef, v::WeakRef) = isequal(w.value, v.value)
==(w::WeakRef, v) = isequal(w.value, v)
==(w, v::WeakRef) = isequal(w, v.value)

function finalizer(o::ANY, f::ANY)
    if isimmutable(o)
        error("objects of type ", typeof(o), " cannot be finalized")
    end
    ccall(:jl_gc_add_finalizer_th, Void, (Ptr{Void}, Any, Any),
          Core.getptls(), o, f)
end
function finalizer{T}(o::T, f::Ptr{Void})
    @_inline_meta
    if isimmutable(T)
        error("objects of type ", T, " cannot be finalized")
    end
    ccall(:jl_gc_add_ptr_finalizer, Void, (Ptr{Void}, Any, Ptr{Void}),
          Core.getptls(), o, f)
end

finalize(o::ANY) = ccall(:jl_finalize_th, Void, (Ptr{Void}, Any,),
                         Core.getptls(), o)

gc(full::Bool=true) = ccall(:jl_gc_collect, Void, (Int32,), full)
gc_enable(on::Bool) = ccall(:jl_gc_enable, Int32, (Int32,), on) != 0

immutable Nullable{T}
    hasvalue::Bool
    value::T

    Nullable() = new(false)
    Nullable(value::T, hasvalue::Bool=true) = new(hasvalue, value)
end

max_stars_repo_name: oxinabox/julia

max_issues_repo_name: oxinabox/julia

max_forks_repo_name: oxinabox/julia