text:
function savepopulation(fname,population)
    fname = fname[end-3:end] == ".csv" ? fname : fname*".csv"
    array = Array{Any,2}(undef,length(population),2+length(population[1].karva))
        for (n,circuit) in enumerate(population)
            array[n,:] = vcat([circuit.karva,circuit.fitness == nothing ? Inf : circuit.fitness],circuit.parameters)
        end
        writedlm(fname,array,';')
end

"""
    loadpopulation(filepath::String)

Load a library of equivalent electrical circuits, to be provided as initial population for the circuitevolution function.

"""

function loadpopulation(filepath)
    population_array = readdlm(filepath,';')
    population_size = size(population_array,1)
    population = Array{Circuit}(undef,population_size)
    cpe_indices = findall(x->typeof(x)==SubString{String},population_array[:,:])
    population_array[cpe_indices[size(population_array,1)+1:end]] = eval.(Meta.parse.(population_array[cpe_indices[population_size+1:end]]))
    for individual in 1:population_size
        population[individual] = Circuit(population_array[individual,1],vcat(population_array[individual,3:end]),population_array[individual,2] == Inf ? nothing : population_array[individual,2])
    end
    return population
end

function readmeasurements(file)
    meansurement_file = readdlm(file,',')
    reals = meansurement_file[:,1]
    imags = meansurement_file[:,2]
    frequencies = meansurement_file[:,3]
    measurements = reals + imags*im
    return measurements,frequencies
end

"""
initialize_circuitlibrary(circuit,parameters)

Initializes a circuit library in the current directory. the circuit
argument should be provided in circuit notation form (e.g. "R1-[C2,R3-[C4,R5]]"),
and the parameters should be provided as an array (e.g. [20,4e-9,3400,4e-6,2500]).
"""

function initialize_circuitlibrary(circuit,parameters=nothing,type = "N",source = "no information provided.")
    @assert !isfile("Circuitlibrary.csv") "A circuitlibrary already exists in the current directory."
    if isnothing(parameters)
        parameters = circuitparameters(circuit)
    end
    circuit_code,circuit_parameters = circuit_to_karva(circuit,parameters)
    circuitlibrary = [Circuit(circuit_code,circuit_parameters,nothing)]
    circuitmetadata = Array{Any}(undef,1,4)
    circuitmetadata[1:4] = [circuit,parameters,type,source]
    savepopulation("Circuitlibrary.csv",circuitlibrary)
    writedlm("circuitmetadata.csv",circuitmetadata,';')
end

"""
add_to_circuitlibrary(circuit,parameters)

Adds a new circuit to an already initialized circuit library. the circuit
argument should be provided in circuit notation form (e.g. "R1-[C2,R3-[C4,R5]]"),
and the parameters should be provided as an array (e.g. [20,4e-9,3400,4e-6,2500]).
"""
function add_to_circuitlibrary(circuit,parameters=nothing,type = "N",source = "no information provided.") #path to circuit library should be able to be provided in function arguments.
    @assert isfile("CircuitLibrary.csv") "A circuitlibrary has not yet been initialized in the current directory."
    if isnothing(parameters)
        parameters = circuitparameters(circuit)
    end
    #convert the circuit to its basic encoding and parameter values.
    circuit_code,circuit_parameters = circuit_to_karva(circuit,parameters)
    circuitobject = Circuit(circuit_code,circuit_parameters,nothing)
    #elongate with random terminals until required length, dictated by the head length.
    #load the circuitlibrary and circuitlibrary metadata files.
    circuitlibrary = loadpopulation("CircuitLibrary.csv")
    metadatalibrary = readdlm("circuitmetadata.csv",';')
    #append the circuit.
    push!(circuitlibrary,circuitobject)
    circuitmetadata = Array{Any}(undef,1,4)
    circuitmetadata[1:4] = [circuit,parameters,type,source]
    metadatalibrary = vcat(metadatalibrary,circuitmetadata)
    #save the circuitlibrary files.
    savepopulation("Circuitlibrary.csv",circuitlibrary)
    writedlm("Circuitmetadata.csv",metadatalibrary,';')
end

function add_to_circuitlibrary_limitedparams(circuit,limitedparams,type = "N",source = "no information provided.") #when only a part of the parameters are known.
    @assert isfile("CircuitLibrary.csv") "A circuitlibrary has not yet been initialized in the current directory."
    parameters = circuitparameters(circuit,limitedparams)
    #convert the circuit to its basic encoding and parameter values.
    circuit_code,circuit_parameters = circuit_to_karva(circuit,parameters)
    circuitobject = Circuit(circuit_code,circuit_parameters,nothing)
    #elongate with random terminals until required length, dictated by the head length.
    #load the circuitlibrary and circuitlibrary metadata files.
    circuitlibrary = loadpopulation("CircuitLibrary.csv")
    metadatalibrary = readdlm("circuitmetadata.csv",';')
    #append the circuit.
    push!(circuitlibrary,circuitobject)
    circuitmetadata = Array{Any}(undef,1,4)
    circuitmetadata[1:4] = [circuit,parameters,type,source]
    metadatalibrary = vcat(metadatalibrary,circuitmetadata)
    #save the circuitlibrary files.
    savepopulation("Circuitlibrary.csv",circuitlibrary)
    writedlm("Circuitmetadata.csv",metadatalibrary,';')
end

function add_encoding_to_circuitlibrary(circuit_code,type = "N",source = "no information provided.")
    circuit_parameters = karva_parameters(circuit_code)
    additional_codings = join(rand("RCLP",17-length(circuit_code))) #assuming a headlength of 8.
    addidional_parameters = karva_parameters(additional_codings)
    karva = circuit_code*additional_codings
    parameters = vcat(circuit_parameters,addidional_parameters)
    circuitlibrary = loadpopulation("CircuitLibrary.csv")
    metadatalibrary = readdlm("circuitmetadata.csv",';')
    circuitobject = Circuit(karva,parameters,nothing)
    push!(circuitlibrary,circuitobject)
    circuitmetadata = Array{Any}(undef,1,4)
    tree =  karva_to_tree(circuit_code,circuit_parameters)
    circuit,parameters_ = tree_to_circuit(tree)
    circuitmetadata[1:4] = [circuit,parameters_,type,source]
    metadatalibrary = vcat(metadatalibrary,circuitmetadata)
     #save the circuitlibrary files.
     savepopulation("Circuitlibrary.csv",circuitlibrary)
     writedlm("Circuitmetadata.csv",metadatalibrary,';')
end

function add_encoding_to_circuitlibrary_with_parameters(circuit_code,parameters,type = "N",source = "no information provided.")
    additional_codings = join(rand("RCLP",17-length(circuit_code))) #assuming a headlength of 8.
    addidional_parameters = karva_parameters(additional_codings)
    karva = circuit_code*additional_codings
    parameters = vcat(parameters,addidional_parameters)
    circuitlibrary = loadpopulation("CircuitLibrary.csv")
    metadatalibrary = readdlm("circuitmetadata.csv",';')
    circuitobject = Circuit(karva,parameters,nothing)
    push!(circuitlibrary,circuitobject)
    circuitmetadata = Array{Any}(undef,1,4)
    tree =  karva_to_tree(circuit_code,circuit_parameters)
    circuit,parameters_ = tree_to_circuit(tree)
    circuitmetadata[1:4] = [circuit,parameters_,type,source]
    metadatalibrary = vcat(metadatalibrary,circuitmetadata)
     #save the circuitlibrary files.
     savepopulation("Circuitlibrary.csv",circuitlibrary)
     writedlm("Circuitmetadata.csv",metadatalibrary,';')
end

function circuitparameters(circuit)
    elements = foldl(replace,["["=>"","]"=>"","-"=>"",","=>""],init = denumber_circuit(circuit))
    return karva_parameters2(elements)
end

function circuitparameters(circuit,limited_params)
    elements = foldl(replace,["["=>"","]"=>"","-"=>"",","=>""],init = denumber_circuit(circuit))
    return karva_parameters3(elements,limited_params)
end

function karva_parameters2(karva)
    parameters = Array{Float64}(undef, length(karva))
    ranges = Dict('R'=>4000,'C'=> 0.0001,'L'=> 1,'+'=> 0,'-'=> 0,'P'=>[4000,1])
    for (e,i) in enumerate(karva)
            if e == 'P'
                Element_values[n] = [ranges[e][1]*rand(),ranges[e][2]*rand()]
            else
            parameters[e] = rand()*ranges[i]
        end
    end
    return parameters
end

function karva_parameters3(karva,limited_params)
    parameters = Array{Any}(undef, length(karva))
    ranges = Dict('R'=>4000,'C'=> 0.0001,'L'=> 1,'+'=> 0,'-'=> 0,'P'=>[4000,1])
    for (e,(i,l)) in enumerate(zip(karva,limited_params))
            if l != 0
                parameters[e] = l
            else
                if i == 'P'
                    Element_values[e] = [ranges[i][1]*rand(),ranges[i][2]*rand()]
                else
                parameters[e] = rand()*ranges[i]
            end
        end
    end
    return parameters
end


text:
# by lazarusA # HIDE
#using GLMakie, Random # HIDE
using CairoMakie, RDatasets, Colors, ColorSchemes
CairoMakie.activate!() #HIDE
let
    airquality = dataset("datasets", "airquality")
    categories = ["Ozone", "Solar.R", "Wind", "Temp"]
    colors = ColorScheme(range(colorant"black", colorant"red",
                length=length(categories)))

    fig = Figure(resolution = (700, 800), font = "sans")
    ax1 = Axis(fig, xlabel = "variable", ylabel = "",
            xticks = (1:length(categories), categories))
    ax2 = Axis(fig, xlabel = "variable", ylabel = "",
            xticks = (1:length(categories), categories))

    for (indx,f) in enumerate(categories)
        datam = []
        for v in airquality[:,f]
            if v !== missing
                push!(datam, v)
            end
        end
        a = fill(indx, length(datam))
        boxplot!(ax1, a, Float32.(datam); whiskerwidth = 1, width = 0.35,
         color = (colors[indx], 0.45), whiskercolor = (colors[indx], 1),
         mediancolor = :black) #show_outliers=false
        violin!(ax2,a,Float32.(datam);width = 0.35,color =(colors[indx], 0.45),
            strokecolor = colors[indx],show_median = true,mediancolor = :black)
    end
    fig[1,1] = ax1
    fig[2,1] = ax2
    fig
    #save("FigAirqualityBox.png", fig, px_per_unit = 2) # HIDE
    save(joinpath(@__DIR__, "output", "BarViolinDataSet.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie", "RDatasets", "Colors", "ColorSchemes"]) # HIDE

ext: jl

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
@defcomp SLRDamages begin
    region = Index()

    y_year = Parameter(index=[time], unit="year")
    y_year_0 = Parameter(unit="year")

    #incoming parameters from SeaLevelRise
    s_sealevel = Parameter(index=[time], unit="m")
    #incoming parameters to calculate consumption per capita after Costs
    cons_percap_consumption = Parameter(index=[time, region], unit="\$/person")
    tct_per_cap_totalcostspercap = Parameter(index=[time,region], unit= "\$/person")
    act_percap_adaptationcosts = Parameter(index=[time, region], unit="\$/person")

    #component parameters
    impmax_maxSLRforadaptpolicySLR = Parameter(index=[region], unit= "m")

    save_savingsrate = Parameter(unit= "%", default=15.00) #pp33 PAGE09 documentation, "savings rate".
    wincf_weightsfactor =Parameter(index=[region], unit="")
    W_SatCalibrationSLR =Parameter(default=1.0) #pp33 PAGE09 documentation, "Sea level impact at calibration sea level rise"
    ipow_SLRIncomeFxnExponent =Parameter(default=-0.30)
    pow_SLRImpactFxnExponent=Parameter(default=0.7333333333333334)
    iben_SLRInitialBenefit=Parameter(default=0.00)
    scal_calibrationSLR = Parameter(default=0.5)
    GDP_per_cap_focus_0_FocusRegionEU = Parameter(unit="\$/person", default=27934.244777382406)

    #component variables
    cons_percap_aftercosts = Variable(index=[time, region], unit = "\$/person")
    gdp_percap_aftercosts = Variable(index=[time, region], unit = "\$/person")

    atl_adjustedtolerablelevelofsealevelrise = Parameter(index=[time,region], unit="m") # meter
    imp_actualreductionSLR = Parameter(index=[time, region], unit="%")
    i_regionalimpactSLR = Variable(index=[time, region], unit="m")

    iref_ImpactatReferenceGDPperCapSLR = Variable(index=[time, region])
    igdp_ImpactatActualGDPperCapSLR=Variable(index=[time, region])

    isatg_impactfxnsaturation = Parameter(unit="unitless")
    isat_ImpactinclSaturationandAdaptationSLR= Variable(index=[time,region])
    isat_per_cap_SLRImpactperCapinclSaturationandAdaptation = Variable(index=[time, region], unit="\$/person")

    rcons_per_cap_SLRRemainConsumption = Variable(index=[time, region], unit = "\$/person") #include?
    rgdp_per_cap_SLRRemainGDP = Variable(index=[time, region], unit = "\$/person")

    function run_timestep(p, v, d, t)

        for r in d.region
            v.cons_percap_aftercosts[t, r] = p.cons_percap_consumption[t, r] - p.tct_per_cap_totalcostspercap[t, r] - p.act_percap_adaptationcosts[t, r]
            v.gdp_percap_aftercosts[t,r]=v.cons_percap_aftercosts[t, r]/(1 - p.save_savingsrate/100)

            if (p.s_sealevel[t]-p.atl_adjustedtolerablelevelofsealevelrise[t,r]) < 0
                v.i_regionalimpactSLR[t,r] = 0
            else
                v.i_regionalimpactSLR[t,r] = p.s_sealevel[t]-p.atl_adjustedtolerablelevelofsealevelrise[t,r]
            end

            v.iref_ImpactatReferenceGDPperCapSLR[t,r]= p.wincf_weightsfactor[r]*((p.W_SatCalibrationSLR + p.iben_SLRInitialBenefit * p.scal_calibrationSLR)*
                (v.i_regionalimpactSLR[t,r]/p.scal_calibrationSLR)^p.pow_SLRImpactFxnExponent - v.i_regionalimpactSLR[t,r] * p.iben_SLRInitialBenefit)

            v.igdp_ImpactatActualGDPperCapSLR[t,r]= v.iref_ImpactatReferenceGDPperCapSLR[t,r]*
                    (v.gdp_percap_aftercosts[t,r]/p.GDP_per_cap_focus_0_FocusRegionEU)^p.ipow_SLRIncomeFxnExponent

            if v.igdp_ImpactatActualGDPperCapSLR[t,r] < p.isatg_impactfxnsaturation
                v.isat_ImpactinclSaturationandAdaptationSLR[t,r] = v.igdp_ImpactatActualGDPperCapSLR[t,r]
            else
                v.isat_ImpactinclSaturationandAdaptationSLR[t,r] = p.isatg_impactfxnsaturation+
                    ((100-p.save_savingsrate)-p.isatg_impactfxnsaturation)*
                    ((v.igdp_ImpactatActualGDPperCapSLR[t,r]-p.isatg_impactfxnsaturation)/
                    (((100-p.save_savingsrate)-p.isatg_impactfxnsaturation)+
                    (v.igdp_ImpactatActualGDPperCapSLR[t,r]- p.isatg_impactfxnsaturation)))
                end
            if v.i_regionalimpactSLR[t,r] < p.impmax_maxSLRforadaptpolicySLR[r]
                v.isat_ImpactinclSaturationandAdaptationSLR[t,r]=v.isat_ImpactinclSaturationandAdaptationSLR[t,r]*(1-p.imp_actualreductionSLR[t,r]/100)
            else
                v.isat_ImpactinclSaturationandAdaptationSLR[t,r]=v.isat_ImpactinclSaturationandAdaptationSLR[t,r]*(1-(p.imp_actualreductionSLR[t,r]/100)* p.impmax_maxSLRforadaptpolicySLR[r] /
                    v.i_regionalimpactSLR[t,r])
                end

                v.isat_per_cap_SLRImpactperCapinclSaturationandAdaptation[t,r] = (v.isat_ImpactinclSaturationandAdaptationSLR[t,r]/100)*v.gdp_percap_aftercosts[t,r]
                v.rcons_per_cap_SLRRemainConsumption[t,r] = v.cons_percap_aftercosts[t,r] - v.isat_per_cap_SLRImpactperCapinclSaturationandAdaptation[t,r]
                v.rgdp_per_cap_SLRRemainGDP[t,r] = v.rcons_per_cap_SLRRemainConsumption[t,r]/(1-p.save_savingsrate/100)

        end

    end
end

# Still need this function in order to set the parameters than depend on
# readpagedata, which takes model as an input. These cannot be set using
# the default keyword arg for now.
function addslrdamages(model::Model)
    SLRDamagescomp = add_comp!(model, SLRDamages)

    SLRDamagescomp[:wincf_weightsfactor] = readpagedata(model, "data/wincf_weightsfactor.csv")
    SLRDamagescomp[:impmax_maxSLRforadaptpolicySLR] = readpagedata(model, "data/impmax_sealevel.csv")

    return SLRDamagescomp
end


text:
using Plots
using StaticArrays
import Base.-, Base.+, Base.*, Base./
using ForwardDiff
using LinearAlgebra

pyplot();

const Scalar = Float64;
const Vector = SVector{2,Float64};
const Label  = Int;

×(u::Vector, v::Vector) = u[1]*v[2] - u[2]*v[1];

const ScalarList = Array{Scalar,1};
const VectorList = Array{Vector,1};
const LabelList  = Array{Label,1};
const LabelListList = Array{LabelList,1};

struct PatchInfo
    name  :: String
    range :: UnitRange{Label}
end;

const PatchInfoList = Array{PatchInfo,1};

struct Mesh
    point      :: VectorList
    centre     :: VectorList
    volume     :: ScalarList
    surface    :: VectorList
    facecentre :: VectorList
    owner      :: LabelList
    neighbor   :: LabelList
    patch      :: PatchInfoList

    cell2point :: LabelListList
    face2point :: LabelListList
end

Mesh() = Mesh(VectorList(), VectorList(), ScalarList(), VectorList(), VectorList(), LabelList(), LabelList(), PatchInfoList(), LabelListList(), LabelListList())

function update!(mesh::Mesh)

    # Recalculate cell centres & volumes
    ncells = size(mesh.cell2point,1)
    empty!(mesh.centre)
    empty!(mesh.volume)
    for c=1:ncells
        pts = mesh.point[mesh.cell2point[c]]

        v = Scalar(0)
        x = Vector(0,0)
        p1 = pts[end]
        for i=1:size(pts,1)-2
            p2 = pts[i]
            p3 = pts[i+1]
            vv = (p2 - p1) × (p3 - p1) / 2
            v += vv
            x += vv*(p1 + p2 + p3)/3
        end
        push!(mesh.volume, v)
        push!(mesh.centre, x/v)
    end

    # Recalculate face sizes & positions
    empty!(mesh.surface)
    empty!(mesh.facecentre)
    nfaces = size(mesh.face2point,1)
    for f=1:nfaces
        pts = mesh.point[mesh.face2point[f]]
        push!(mesh.surface, Vector(pts[2][2]-pts[1][2], pts[1][1]-pts[2][1]))
        push!(mesh.facecentre, (pts[1]+pts[2])/2)
    end
end

function cartesian_mesh(nx, ny)
    Δx, Δy = 1.0/nx, 1.0/ny

    pid(i,j) = i + j*(nx+1) + 1
    cid(i,j) = i + (j-1)*nx

    mesh = Mesh()

    # Mesh points
    for j=0:ny, i=0:nx
        push!(mesh.point, Vector(i*Δx, j*Δy))
    end

    # Mesh cells
    for j=1:ny, i=1:nx
        pts = [pid(i-1,j-1), pid(i,j-1), pid(i,j), pid(i-1,j)] |> LabelList
        push!(mesh.cell2point, pts)
    end

    # Internal faces
    for j=1:ny, i=1:nx-1
        push!(mesh.owner, cid(i,j))
        push!(mesh.neighbor, cid(i+1,j))
        push!(mesh.face2point, [pid(i,j-1), pid(i,j)])
    end

    for j=1:ny-1, i=1:nx
        push!(mesh.owner, cid(i,j))
        push!(mesh.neighbor, cid(i,j+1))
        push!(mesh.face2point, [pid(i,j), pid(i-1,j)])
    end

    # Boundary patches
    # - bottom
    j=0
    start = size(mesh.owner,1)+1
    for i=1:nx
        push!(mesh.owner, cid(i,j+1))
        push!(mesh.face2point, [pid(i-1,j), pid(i,j)])
    end
    push!(mesh.patch, PatchInfo("bottom", UnitRange(start, size(mesh.owner,1))))

    # - right
    i=nx
    start = size(mesh.owner,1)+1
    for j=1:ny
        push!(mesh.owner, cid(i,j))
        push!(mesh.face2point, [pid(i,j-1), pid(i,j)])
    end
    push!(mesh.patch, PatchInfo("right", UnitRange(start, size(mesh.owner,1))))

    # - top
    j=ny
    start = size(mesh.owner,1)+1
    for i=1:nx
        push!(mesh.owner, cid(i,j))
        push!(mesh.face2point, [pid(i,j), pid(i-1,j)])
    end
    push!(mesh.patch, PatchInfo("top", UnitRange(start, size(mesh.owner,1))))

    # - left
    i=0
    start = size(mesh.owner,1)+1
    for j=1:ny
        push!(mesh.owner, cid(i+1,j))
        push!(mesh.face2point, [pid(i,j), pid(i,j-1)])
    end
    push!(mesh.patch, PatchInfo("left", UnitRange(start, size(mesh.owner,1))))

    update!(mesh)

    return mesh
end

internal_faces(m::Mesh) = UnitRange(1,size(m.neighbor,1))

boundary_patches(m::Mesh) = UnitRange(1,size(m.patch,1))

patch_names(m::Mesh) = [ m.patch[p].name for p in boundary_patches(m)]

patch_by_name(m::Mesh, name::String) = findfirst(x->x==name, patch_names(m))

patch_faces(m::Mesh, patch::Label) = m.patch[patch].range

m2 = cartesian_mesh(2,2);

for p in m2.point
    println("p = ", p)
end

for c in m2.centre
    println("c = ", c)
end

print("volume[:] = ", m2.volume)

for f in internal_faces(m2)
    println("S=", m2.surface[f], "\t o=", m2.owner[f], "\t n=",m2.neighbor[f])
end

for p in boundary_patches(m2)
    println("PATCH ", m2.patch[p].name)
    for f in patch_faces(m2, p)
        println("S=", m2.surface[f], "\t o=", m2.owner[f])
    end
    println()
end

mutable struct Field
    values
    mesh
    boundaries
end

Field(m::Mesh) = Field( zeros(length(m.centre)), m, Dict{String,Any}());
Field(values, m::Mesh) = Field(values, m, Dict{String,Any}());

for op ∈ [:+, :-, :*, :/]
    @eval begin
        function ($op)(a::Field, b::Field)
            @assert a.mesh == b.mesh
            Field(broadcast($op,a.values,b.values), a.mesh)
        end
    end
end

a2 = Field(m2); a2.values = ones(4);
b2 = Field(m2); b2.values = 3*ones(4);

c2 = a2 * b2;
c2.values

for op ∈ [:+, :-, :*, :/]
    @eval begin
        function ($op)(a::Field, b::Number)
            Field(broadcast($op,a.values,b), a.mesh)
        end
    end
end

for op ∈ [:+, :-, :*, :/]
    @eval begin
        function ($op)(a::Field, b::Array{Float64,1})
            @assert size(a.values) == size(b)
            Field(broadcast($op,a.values,b), a.mesh)
        end
    end
end

-(a::Field) = Field(-a.values, a.mesh);

c2 = a2 / 3;
c2.values

function ←(a::Field, b::Field)
    @assert a.mesh == b.mesh
    a.values = copy(b.values)
end

b2.values = 3*ones(4)
c2 ← b2;
d2 = b2;
b2.values[1] = 5;
println(c2.values)
println(d2.values)

struct BoundaryCondition
    c1    :: ScalarList
    c2    :: ScalarList
    field :: Field
    patch :: Label
end

function boundary_coeffs(bc::BoundaryCondition, f::Label)
    field = bc.field
    mesh = field.mesh
    start = first(mesh.patch[bc.patch].range)
    i = f + 1 - start

    c1 = bc.c1[i]
    c2 = bc.c2[i]

    return c1, c2
end

function boundary_value(bc::BoundaryCondition, f::Label)
    field = bc.field
    mesh = field.mesh
    u1 = field.values[mesh.owner[f]]
    c1, c2 = boundary_coeffs(bc, f)
    return c1 + c2*u1
end

function set_dirichlet_patch!(field::Field, name::String, val::Scalar)
    mesh = field.mesh
    p = patch_by_name(mesh, name)
    c1 = [val for f in patch_faces(mesh, p)]
    c2 = zero(c1)
    field.boundaries[name] = BoundaryCondition(c1, c2, field, p)
end;

function set_neumann_patch!(field::Field, name::String)
    mesh = field.mesh
    p = patch_by_name(mesh, name)
    c2 = [1.0 for f in patch_faces(mesh, p)]
    c1 = zero(c2)
    field.boundaries[name] = BoundaryCondition(c1, c2, field, p)
end;

m3 = cartesian_mesh(3,3);

u3 = Field(m3);
set_dirichlet_patch!(u3, "left", 1.0);
set_dirichlet_patch!(u3, "bottom", 0.0);
set_neumann_patch!(u3, "right");
set_neumann_patch!(u3, "top");

v = Vector(1, 0.5);

# Tok typu upwind
H(ul, ur, s) = dot(v,s)>0 ? dot(v,s)*ul : dot(v,s)*ur

function R(u::Field)
    mesh = u.mesh
    r = zero(u.values)

    for f in internal_faces(mesh)
        o, n = mesh.owner[f], mesh.neighbor[f]
        S = mesh.surface[f]
        flux = H(u.values[o], u.values[n], S)
        r[o] += flux
        r[n] -= flux
    end

    for p in boundary_patches(mesh)
        name = mesh.patch[p].name
        bc = u.boundaries[name]

        for f in patch_faces(mesh, p)
            o = mesh.owner[f]
            S = mesh.surface[f]
            ub = boundary_value(bc, f)
            flux = H(u.values[o], ub, S)
            r[o] += flux
        end
    end

    r ./= mesh.volume

    return r
end

R(u3)

m = cartesian_mesh(50,50);

u = Field(m)
set_dirichlet_patch!(u, "left", 1.0);
set_dirichlet_patch!(u, "bottom", 0.0);
set_neumann_patch!(u, "right");
set_neumann_patch!(u, "top");

Δt = 0.01
for t = 0 : Δt : 2
    u ← u - Δt * R(u)
end

plot(reshape(u.values, 50,50)', st=:contour, color=:lightrainbow, fill=true, aspect_ratio=:equal)

using SparseArrays;

struct Equation
    A
    x
    b
end

# Scitani a odcitani rovnic
for op ∈ [:+, :-]
    @eval begin
        function ($op)(e1::Equation, e2::Equation)
            @assert e1.x == e2.x
            Equation(($op)(e1.A, e2.A), e1.x, $(op)(e1.b, e2.b))
        end
    end
end

# Pricteni/odecteni pole o spravne velikosti
for op ∈ [:+, :-]
    @eval begin
        function ($op)(e1::Equation, f::Array{Float64,1})
            @assert size(e1.x) == size(f)
            Equation(copy(e1.A), e1.x, $(op)(e1.b, f))
        end
    end
end

# Vynasobeni rovnice konstantou
*(a::Float64, eq::Equation) = Equation(a*eq.A, eq.x, a*eq.b);

function solve!(eqn::Equation)
    x = eqn.A \ eqn.b
    for i in eachindex(eqn.x)
        eqn.x[i] = -x[i]
    end
end;

function ddt(f::Field,Δt)
    n = length(f.values)
    Equation(I/Δt, f.values, -f.values/Δt)
end;

ddt(u3,0.1)

u3.values = zeros(9)
for t=1:10
    solve!(ddt(u3,0.1)-ones(9))
    print(u3.values[1], "\t")
end

function ∇(v::Vector, u::Field)
    A = spzeros(length(u.values),length(u.values))
    b = zero(u.values)

    msh = u.mesh

    for f in internal_faces(msh)
        o = msh.owner[f]
        n = msh.neighbor[f]
        S = msh.surface[f]
        uo = u.values[o]
        un = u.values[n]

        ϕ = dot(v,S)
        α = max(ϕ, 0.0)
        β = min(ϕ, 0.0)

        A[o,o] += α / msh.volume[o]
        A[o,n] += β / msh.volume[o]

        A[n,o] -= α / msh.volume[n]
        A[n,n] -= β / msh.volume[n]
    end

    for p in boundary_patches(msh)
        name = msh.patch[p].name
        bc = u.boundaries[name]

        for f in patch_faces(msh, p)
            o = msh.owner[f]
            S = msh.surface[f]

            ϕ = dot(v,S)
            α = max(ϕ, 0.0)
            β = min(ϕ, 0.0)

            # ub = c1 + c2*u[o], flux=alpha*u[o] + beta*ub
            c1, c2 = boundary_coeffs(bc, f)

            A[o,o] += (α + β*c2) / msh.volume[o]
            b[o]   += β*c1 / msh.volume[o]
        end
    end

    return Equation(A, u.values, b)
end

Matrix(∇(v,u3).A)

u.values = zero(u.values)

Δt = 0.05
for t=0:Δt:0.5
    solve!( ddt(u,Δt) + ∇(v,u) );
end

plot(reshape(u.values, 50,50)', st=:contour, color=:lightrainbow, fill=true, aspect_ratio=:equal)

function Δ(μ::Float64, u::Field)
    mesh = u.mesh

    A = spzeros(length(u.values),length(u.values))
    b = zero(u.values)

    for f in internal_faces(mesh)
        o, n = mesh.owner[f], mesh.neighbor[f]
        xo, xn = mesh.centre[o], mesh.centre[n]
        S = mesh.surface[f]

        g = μ * norm(S) / norm(xn - xo)

        A[o,o] -= g / mesh.volume[o]
        A[o,n] += g / mesh.volume[o]

        A[n,o] += g / mesh.volume[n]
        A[n,n] -= g / mesh.volume[n]
    end

    for p in boundary_patches(mesh)
        name = mesh.patch[p].name
        bc = u.boundaries[name]

        for f in patch_faces(mesh, p)
            o = mesh.owner[f]
            S = mesh.surface[f]
            xo, xf = mesh.centre[o], mesh.facecentre[f]
            c1, c2 = boundary_coeffs(bc, f)

            # normal derivative at boundary dudn = (ub - uin) / delta
            # with ub = c[1] + c[2]*uin, i.e.
            # dudn = c[1]/delta + (c[2]-1)/delta*uin

            δ = dot(xf - xo, S)/norm(S)
            A[o,o] += μ * (c2-1)/δ * norm(S) / mesh.volume[o]
            b[o] += μ * c1/δ * norm(S) / mesh.volume[o]
        end
    end

    Equation(A, u.values, b)
end

D=Matrix(Δ(1.0,u3).A)

A = Matrix((ddt(u3,Δt) + ∇(v,u3) - Δ(1.e-1,u3)).A)

u.values = zero(u.values)

Δt = 0.1
for t=0:Δt:0.5
    solve!( ddt(u,Δt) + ∇(v,u) - Δ(1.e-1,u) );
end

plot(reshape(u.values, 50,50)', st=:contour, color=:lightrainbow, fill=true, aspect_ratio=:equal)

max_stars_repo_name: furstj/NSP

max_issues_repo_name: furstj/NSP

max_forks_repo_name: furstj/NSP


text:
# Autogenerated wrapper script for MariaDB_Connector_ODBC_jll for x86_64-linux-musl
export libmaodbc

using iODBC_jll
using Libiconv_jll
using unixODBC_jll
using OpenSSL_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("MariaDB_Connector_ODBC")
JLLWrappers.@declare_library_product(libmaodbc, "libmaodbc.so")
function __init__()
    JLLWrappers.@generate_init_header(iODBC_jll, Libiconv_jll, unixODBC_jll, OpenSSL_jll, Zlib_jll)
    JLLWrappers.@init_library_product(
        libmaodbc,
        "lib/mariadb/libmaodbc.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Silverman's rule of thumb for KDE bandwidth selection
function bandwidth(data::Vector, alpha::Float64 = 0.9)
    # Determine length of data
    ndata = length(data)

    # Calculate width using variance and IQR
    var_width = std(data)
    q25, q75 = quantile(data, [0.25, 0.75])
    quantile_width = (q75 - q25) / 1.34

    # Deal with edge cases with 0 IQR or variance
    width = min(var_width, quantile_width)
    if width == 0.0
        if var_width == 0.0
            width = 1.0
        else
            width = var_width
        end
    end

    # Set bandwidth using Silverman's rule of thumb
    return alpha * width * ndata^(-0.2)
end

# Store both grid and density for KDE over the real line
immutable UnivariateKDE
    x::Vector{Float64}
    density::Vector{Float64}
end

# Store both grid and density for KDE over R2
immutable BivariateKDE
    x::Vector{Float64}
    y::Vector{Float64}
    density::Matrix{Float64}
end

# Algorithm AS 176 for calculating univariate KDE
function kde{T <: Real}(data::Vector{T},
                        window::Real = bandwidth(data),
                        npoints::Integer = 2048)
    # Determine length of data
    ndata = length(data)

    # Check that the window is a positive constant
    if window <= 0.0
        error("Window must be positive")
    end

    # Find interval that will contain almost all mass
    dlo = min(data) - 3 * window
    dhi = max(data) + 3 * window

    # Check that interval for estimation is valid
    if dlo >= dhi
        error("Interval must be specified as lower bound, upper bound")
    end

    # Set up a grid for discretized data
    grid = zeros(Float64, npoints)

    # Define some more constants
    step = (dhi - dlo) / npoints
    ainc = 1.0 / (ndata * step)
    npoints2 = fld(npoints, 2)
    hw = window / step
    fac1 = 32.0 * (atan(1.0) * hw / npoints)^2

    # Discretize the data using a histogram
    dlo1 = dlo - step
    for i in 1:ndata
        j = iround(fld(data[i] - dlo1, step))
        if j >= 1 && j <= npoints
            grid[j] = grid[j] + ainc
        end
    end

    # Transform to Fourier basis
    ft = rfft(grid)

    # Find transform of KDE by convolving grid with the
    # Fourier transform of a Gaussian kernel
    for j = 2:length(ft)
        ft[j] *= exp(-fac1 * (j-1)^2)
    end

    # Invert the Fourier transform to get the KDE
    density = irfft(ft, npoints)

    # Fix any noise that crept in
    for j in 1:npoints
        if density[j] < 0.0
            density[j] = 0.0
        end
    end

    # Expand the grid over which KDE was calculated
    x = [dlo:step:(dhi - step)]

    return UnivariateKDE(x, density)
end

# Algorithm from MASS Chapter 5 for calculating 2D KDE
function kde{S <: Real, T <: Real}(x::Vector{S},
                                   y::Vector{T};
                                   h::Real = NaN,
                                   resolution::Integer = 25)
    n = length(x)

    if length(y) != n
        error("x and y must have the same length")
    end

    if isnan(h)
        h1 = bandwidth(x)
        h2 = bandwidth(y)
    else
        h1 = h
        h2 = h
    end

    min_x, max_x = minmax(x)
    min_y, max_y = minmax(y)

    grid_x = [min_x:((max_x - min_x) / (resolution - 1)):max_x]
    grid_y = [min_y:((max_y - min_y) / (resolution - 1)):max_y]

    mx = Array(Float64, resolution, n)
    my = Array(Float64, resolution, n)
    for i in 1:resolution
        for j in 1:n
            mx[i, j] = pdf(Normal(), (grid_x[i] - x[j]) / h1)
            my[i, j] = pdf(Normal(), (grid_y[i] - y[j]) / h2)
        end
    end

    z = A_mul_Bt(mx, my)
    for i in 1:(resolution^2)
        z[i] /= (n * h1 * h2)
    end

    return BivariateKDE(grid_x, grid_y, z)
end


text:
function plane_transformation_strain(theta)
    return [theta / 3.0, theta / 3.0, 0.0]
end

function plane_strain_transformation_stress(lambda, mu, theta)
    pt = (lambda + 2mu / 3) * theta
    transfstress = [pt, pt, 0.0]
    return transfstress
end

function bulk_transformation_rhs(basis, quad, transfstress, jac)
    dim = dimension(basis)
    nf = number_of_basis_functions(basis)
    ndofs = dim * nf
    rhs = zeros(ndofs)
    vectosymmconverter = vector_to_symmetric_matrix_converter()
    detjac = prod(jac)

    for (p, w) in quad
        grad = transform_gradient(gradient(basis, p), jac)
        NK = sum([make_row_matrix(vectosymmconverter[k], grad[:, k]) for k = 1:dim])
        rhs .+= NK' * transfstress * detjac * w
    end
    return rhs
end

function interface_transformation_rhs(basis, quad, normals, transfstress, cellmap)
    numqp = length(quad)
    @assert size(normals)[2] == numqp
    dim = dimension(basis)
    nf = number_of_basis_functions(basis)
    ndofs = dim * nf
    rhs = zeros(ndofs)

    vectosymmconverter = vector_to_symmetric_matrix_converter()
    scalearea = scale_area(cellmap, normals)

    for qpidx = 1:numqp
        p, w = quad[qpidx]

        vals = basis(p)
        NI = interpolation_matrix(vals, dim)

        normal = normals[:, qpidx]
        NK = sum([normal[k] * vectosymmconverter[k]' for k = 1:dim])

        rhs .+= NI' * NK * transfstress * scalearea[qpidx] * w
    end
    return rhs
end

function interface_transformation_component_rhs(
    basis,
    quad,
    components,
    normals,
    transfstress,
    cellmap,
)

    numqp = length(quad)
    @assert size(normals)[2] == size(components)[2] == numqp
    dim = dimension(basis)
    nf = number_of_basis_functions(basis)
    ndofs = dim * nf
    rhs = zeros(ndofs)

    vectosymmconverter = vector_to_symmetric_matrix_converter()
    scalearea = scale_area(cellmap, normals)

    for qpidx = 1:numqp
        p, w = quad[qpidx]

        vals = basis(p)
        NI = interpolation_matrix(vals, dim)

        normal = normals[:, qpidx]
        NK = sum([normal[k] * vectosymmconverter[k]' for k = 1:dim])

        component = components[:,qpidx]
        projector = component*component'

        rhs .+= NI' * projector * NK * transfstress * scalearea[qpidx] * w
    end
    return rhs
end

function face_traction_transformation_rhs(basis, quad, normal, transfstress, facescale)
    dim = dimension(basis)
    @assert length(normal) == dim
    nf = number_of_basis_functions(basis)
    ndofs = dim * nf
    rhs = zeros(ndofs)

    vectosymmconverter = vector_to_symmetric_matrix_converter()
    NK = sum([normal[k] * vectosymmconverter[k]' for k = 1:dim])

    for (p, w) in quad
        vals = basis(p)
        NI = interpolation_matrix(vals, dim)

        rhs .+= NI' * NK * transfstress * facescale * w
    end
    return rhs
end

function face_traction_component_transformation_rhs(
    basis,
    quad,
    component,
    normal,
    transfstress,
    facescale,
)

    dim = dimension(basis)
    @assert length(normal) == dim
    nf = number_of_basis_functions(basis)
    ndofs = dim * nf
    rhs = zeros(ndofs)

    projector = component * component'
    vectosymmconverter = vector_to_symmetric_matrix_converter()
    NK = sum([normal[k] * vectosymmconverter[k]' for k = 1:dim])

    for (p, w) in quad
        vals = basis(p)
        NI = interpolation_matrix(vals, dim)

        rhs .+= NI' * projector * NK * transfstress * facescale * w
    end
    return rhs
end


text:
import MathProgBase
export solve!

# semantics:
# * we *load* problem data from Convex into the MathProgBase model
# * we *populate* objects in Convex with their optimal values

function solve!(problem::Problem,
                s::MathProgBase.AbstractMathProgSolver;
                kwargs...)
    # TODO: warn about wiping out old model if warmstart=true?
    problem.model = MathProgBase.ConicModel(s)
    return solve!(problem; kwargs...)
end

function solve!(problem::Problem;
                warmstart=false,
                check_vexity=true,
                verbose=true)

    if problem.model === nothing
        throw(ArgumentError(
            "The provided problem hasn't been initialized with a conic model.
            You can resolve this by passing in a `AbstractMathProgSolver`. For example,
            ```
            using ECOS
            solve!(problem, ECOSSolver())
            ```"
        ))
    end

    if check_vexity
        vex = vexity(problem)
    end

    c, A, b, cones, var_to_ranges, vartypes, conic_constraints, id_to_variables, conic_constr_to_constr = conic_problem(problem)

    # load MPB conic problem
    m = problem.model
    load_problem!(m, c, A, b, cones, vartypes)
    if warmstart
        set_warmstart!(m, problem, length(c), var_to_ranges, id_to_variables)
    end
    # optimize MPB conic problem
    MathProgBase.optimize!(m)

    # populate the status, the primal (and possibly dual) solution
    # and the primal (and possibly dual) variables with values
    populate_solution!(m, problem, var_to_ranges, conic_constraints, id_to_variables, conic_constr_to_constr)
    if problem.status != :Optimal && verbose
        @warn "Problem status $(problem.status); solution may be inaccurate."
    end
end

function set_warmstart!(m::MathProgBase.AbstractConicModel,
                        problem::Problem,
                        n::Int, # length of primal (conic) solution
                        var_to_ranges,
                        id_to_variables)
     # use previously cached solution, if any,
     try
         primal = problem.solution.primal
     catch
         @warn "Unable to use cached solution to warmstart problem.
                (Perhaps this is the first time you're solving this problem?)
                Warmstart may be ineffective."
         primal = zeros(n)
     end
     if length(primal) != n
         @warn "Unable to use cached solution to warmstart problem.
                (Perhaps the number of variables or constraints in the problem have changed since you last solved it?)
                Warmstart may be ineffective."
         primal = zeros(n)
     end

     # grab any variables whose values the user may be trying to set
     load_primal_solution!(primal, var_to_ranges, id_to_variables)

     # notify the model that we're trying to warmstart
     try
         MathProgBase.setwarmstart!(m, primal)
     catch
         @warn "Unable to warmstart solution.
               (Perhaps the solver doesn't support warm starts?)
               Using a cold start instead."
     end
     m
end

function load_problem!(m::MathProgBase.AbstractConicModel, c, A, b, cones, vartypes)
    # no conic constraints on variables
    var_cones = fill((:Free, 1:size(A, 2)),1)
    MathProgBase.loadproblem!(m, vec(Array(c)), A, vec(Array(b)), cones, var_cones)

    # add integer and binary constraints on variables
    if !all(==(:Cont), vartypes)
        try
            MathProgBase.setvartype!(m, vartypes)
        catch
            error("model $(typeof(m)) does not support variables of some of the following types: $(unique(vartypes))")
        end
    end
    m
end

function populate_solution!(m::MathProgBase.AbstractConicModel,
                            problem::Problem,
                            var_to_ranges,
                            conic_constraints,
                            id_to_variables,
                            conic_constr_to_constr)
    dual = try
        MathProgBase.getdual(m)
    catch
        fill(NaN, MathProgBase.numconstr(m))
    end

    solution = try
        MathProgBase.getsolution(m)
    catch
        fill(NaN, MathProgBase.numvar(m))
    end

    objective = try
        MathProgBase.getobjval(m)
    catch
        NaN
    end

    if any(isnan, dual)
        problem.solution = Solution(solution, MathProgBase.status(m), objective)
    else
        problem.solution = Solution(solution, dual, MathProgBase.status(m), objective)
    end

    populate_variables!(problem, var_to_ranges, id_to_variables)

    if problem.solution.has_dual
        populate_duals!(conic_constraints, problem.solution.dual, conic_constr_to_constr)
    end

    # minimize -> maximize
    if (problem.head == :maximize)
        problem.solution.optval = -problem.solution.optval
    end

    # Populate the problem with the solution
    problem.optval = problem.solution.optval
    problem.status = problem.solution.status

    problem
end

function populate_variables!(problem::Problem, var_to_ranges::Dict{UInt64, Tuple{Int, Int}}, id_to_variables)
    x = problem.solution.primal
    for (id, (start_index, end_index)) in var_to_ranges
        var = id_to_variables[id]
        sz = var.size
        if var.sign != ComplexSign()
            var.value = reshape(x[start_index:end_index], sz[1], sz[2])
            if sz == (1, 1)
                var.value = var.value[1]
            end
        else
            real_value = reshape(x[start_index:start_index + div(end_index-start_index+1,2)-1], sz[1], sz[2])
            imag_value = reshape(x[start_index + div(end_index-start_index+1,2):end_index], sz[1], sz[2])
            var.value = real_value + im*imag_value
            if sz == (1, 1)
                var.value = var.value[1]
            end
        end
    end
end

# populates the solution vector from the .value fields of variables
# for use in warmstarting
# TODO: it would be super cool to grab the other expressions that appear in the primal solution vector,
# get their `expression_to_range`,
# and populate them too using `evaluate`
function load_primal_solution!(primal::Array{Float64,1}, var_to_ranges::Dict{UInt64, Tuple{Int, Int}}, id_to_variables)
    for (id, (start_index, end_index)) in var_to_ranges
        var = id_to_variables[id]
        if var.value !== nothing
            sz = size(var.value)
            if length(sz) <= 1
                primal[start_index:end_index] = var.value
            else
                primal[start_index:end_index] = reshape(var.value, sz[1]*sz[2], 1)
            end
        end
    end
end

function populate_duals!(constraints::Array{ConicConstr}, dual::Vector, conic_constr_to_constr)
    constr_index = 1
    for constraint in constraints
        # conic_constr_to_constr only has keys for conic constraints with a single objective
        # so this will work
        if haskey(conic_constr_to_constr, constraint)
            sz = constraint.sizes[1]
            c = conic_constr_to_constr[constraint]
            c.dual = reshape(dual[constr_index:constr_index+sz-1], c.size)
            if c.size == (1, 1)
                c.dual = c.dual[1]
            end
            constr_index += sz
        else
            for i = 1:length(constraint.objs)
                sz = constraint.sizes[i]
                constr_index += sz
            end
        end
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct FileStatuses <: SwaggerModel
    FileStatus::Any # spec type: Union{ Nothing, Vector{FileStatusProperties} } # spec name: FileStatus

    function FileStatuses(;FileStatus=nothing)
        o = new()
        validate_property(FileStatuses, Symbol("FileStatus"), FileStatus)
        setfield!(o, Symbol("FileStatus"), FileStatus)
        o
    end
end # type FileStatuses

const _property_map_FileStatuses = Dict{Symbol,Symbol}(Symbol("FileStatus")=>Symbol("FileStatus"))
const _property_types_FileStatuses = Dict{Symbol,String}(Symbol("FileStatus")=>"Vector{FileStatusProperties}")
Base.propertynames(::Type{ FileStatuses }) = collect(keys(_property_map_FileStatuses))
Swagger.property_type(::Type{ FileStatuses }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_FileStatuses[name]))}
Swagger.field_name(::Type{ FileStatuses }, property_name::Symbol) =  _property_map_FileStatuses[property_name]

function check_required(o::FileStatuses)
    true
end

function validate_property(::Type{ FileStatuses }, name::Symbol, val)
end


text:
## quadrature.jl : nystrom integration methods

abstract type QuadratureRule end

struct GaussLegendre <: QuadratureRule end
struct EOLE <: QuadratureRule end
struct Simpson <: QuadratureRule end
struct Midpoint <: QuadratureRule end
struct Trapezoidal <: QuadratureRule end

# get nodes and weights of Gauss-Legendre quadrature on [a,b]
function get_nodes_and_weights(n,a,b,q::GaussLegendre)
    nodes, weights = gausslegendre(n)
    weights = (b-a)/2*weights
    nodes = (b-a)/2*nodes .+ (a+b)/2
    return nodes, weights
end

# get nodes and weights of structured grid on [a,b]
function get_nodes_and_weights(n,a,b,q::EOLE)
    nodes = range(a,stop = b,length = n)
    weights = (b-a)/n*fill(1.0, size(nodes))
    return nodes, weights
end

# get nodes and weights of Simpson's rule on [a,b]
function get_nodes_and_weights(n,a,b,q::Simpson)
    n % 2 == 0 || begin
        warn("to use Simpson's rule, n must be even (received $(n)). I will continue with n = $(n+1)")
        n += 1
    end
    Δx = (b-a)/n
    nodes = a:Δx:b
    weights = repeat(2:2:4, outer=Int(n/2))
    weights[1] = 1
    push!(weights,1)
    weights *= Δx/3
    return nodes, weights
end

# get nodes and weights of Midpoint rule on [a,b]
function get_nodes_and_weights(n,a,b,q::Midpoint)
    Δx = (b-a)/n
    nodes = a+Δx/2:Δx:b-Δx/2
    weights = Δx*fill(1.0, size(nodes))
    return nodes, weights
end

# get nodes and weights of Trapezoidal rule on [a,b]
function get_nodes_and_weights(n,a,b,q::Trapezoidal)
    Δx = (b-a)/n
    nodes = a:Δx:b
    weights = Δx*fill(1.0, size(nodes))
    weights[1] /= 2
    weights[end] /= 2
    return nodes, weights
end


text:
export  Measure,
        ProductMeasure,
        LegendreMeasure,
        JacobiMeasure,
        LaguerreMeasure,
        genLaguerreMeasure,
        HermiteMeasure,
        genHermiteMeasure,
        GaussMeasure,
        Beta01Measure,
        GammaMeasure,
        Uniform01Measure,
        Uniform_11Measure,
        LogisticMeasure,
        MeixnerPollaczekMeasure

Type_for_domain = Tuple{Float64, Float64}
Type_for_function = Tuple{<:Real, <:Real}

struct Measure <: AbstractMeasure
    name::String
    w::Function
    dom::Type_for_domain
    symmetric::Bool
    pars::Dict
    function Measure(name::String, w::Function, dom::Type_for_function, symm::Bool, d::Dict=Dict())
        !(dom[1] < dom[2]) && throw(DomainError(dom, "invalid domain bounds specified"))
        new(lowercase(name),w, dom, symm, d)
    end
end

struct ProductMeasure <: AbstractMeasure
    w::Function
    measures::Vector{<:AbstractMeasure}
end

# constructor for classic distributions
struct LegendreMeasure <: AbstractCanonicalMeasure
    w::Function
    dom::Type_for_domain
    symmetric::Bool
    function LegendreMeasure()
        new(w_legendre,(-1.,1.),true)
    end
end

struct JacobiMeasure <: AbstractCanonicalMeasure
    w::Function
    dom::Type_for_domain
    symmetric::Bool
    ashapeParameter::Float64
    bshapeParameter::Float64

    function JacobiMeasure(shape_a::Real, shape_b::Real)
        any
        shape_a <= -1 && throw(DomainError(shape_a, "shape parameter a must be > -1"))
        shape_b <= -1 && throw(DomainError(shape_b, "shape parameter b must be > -1"))
        new(build_w_jacobi(shape_a,shape_b), (-1,1), isapprox(shape_a, shape_b), shape_a, shape_b)
    end
end

struct LaguerreMeasure <: AbstractCanonicalMeasure
    w::Function
    dom::Type_for_domain
    symmetric::Bool

    function LaguerreMeasure()
        new(w_laguerre, (0,Inf), false)
    end
end

struct genLaguerreMeasure <: AbstractCanonicalMeasure
    w::Function
    dom::Type_for_domain
    symmetric::Bool
    shapeParameter::Float64

    function genLaguerreMeasure(shape::Real)
        shape <= -1 && throw(DomainError(shape, "invalid shape parameter"))
        new(build_w_genlaguerre(shape), (0,Inf), false, shape)
    end
end

struct HermiteMeasure <: AbstractCanonicalMeasure
    w::Function
    dom::Type_for_domain
    symmetric::Bool

    function HermiteMeasure()
        new(w_hermite, (-Inf,Inf), true)
    end
end

struct genHermiteMeasure <: AbstractCanonicalMeasure
    w::Function
    dom::Type_for_domain
    symmetric::Bool
    muParameter::Float64

    function genHermiteMeasure(mu::Real)
        mu <= -0.5 && throw(DomainError(mu, "invalid parameter value (must be > - 0.5)"))
        new(build_w_genhermite(mu), (-Inf,Inf), true, mu)
    end
end

struct MeixnerPollaczekMeasure <: AbstractCanonicalMeasure
    w::Function
    dom::Type_for_domain
    symmetric::Bool
    λParameter::Float64
    ϕParameter::Float64

    function MeixnerPollaczekMeasure(λ::Real, ϕ::Real)
        λ <= 0 && throw(DomainError(λ, "λ has to be positive"))
        !(0 < ϕ < pi) && throw(DomainError(ϕ, "ϕ has to be between 0 and pi"))
        new(build_w_meixner_pollaczek(λ, ϕ), (-Inf,Inf), false, λ, ϕ)
    end
end

struct GaussMeasure <: AbstractCanonicalMeasure
    w::Function
    dom::Type_for_domain
    symmetric::Bool

    function GaussMeasure()
        new(w_gaussian, (-Inf, Inf), true)
    end
end

struct Uniform01Measure <: AbstractCanonicalMeasure
    w::Function
    dom::Type_for_domain
    symmetric::Bool

    function Uniform01Measure()
        new(w_uniform01, (0,1), true)
    end
end

struct Uniform_11Measure <: AbstractCanonicalMeasure
    w::Function
    dom::Type_for_domain
    symmetric::Bool

    function Uniform_11Measure()
        new(w_uniform_11, (-1,1), true)
    end
end

struct Beta01Measure <: AbstractCanonicalMeasure
    w::Function
    dom::Type_for_domain
    symmetric::Bool
    ashapeParameter::Float64
    bshapeParameter::Float64

    function Beta01Measure(a, b)
        a <= 0 && throw(DomainError(a, "shape parameter a must be positive"))
        b <= 0 && throw(DomainError(b, "shape parameter b must be positive"))
        new(build_w_beta(a,b), (0,1), isapprox(a,b), a, b)
    end
end

struct GammaMeasure <: AbstractCanonicalMeasure
    w::Function
    dom::Type_for_domain
    symmetric::Bool
    shapeParameter::Float64
    rateParameter::Float64

    function GammaMeasure(shape::Real, rate::Real)
        shape <= 0 && throw(DomainError(shape, "shape parameter needs to be positive"))
        rate != 1 && throw(DomainError(rate, "rate must be unity (currently!)"))
        new(build_w_gamma(shape), (0,Inf), false, shape, rate)
    end
end

struct LogisticMeasure <: AbstractCanonicalMeasure
    w::Function
    dom::Type_for_domain
    symmetric::Bool

    function LogisticMeasure()
        new(w_logistic, (-Inf,Inf), true)
    end
end


text:
#using Gallery
#using GalleryWaveguide

#OBS: Only needed to run the debug
push!(LOAD_PATH, string(@__DIR__, "/../src/gallery_extra"))
using GalleryWaveguide
push!(LOAD_PATH, string(@__DIR__, "/../src/gallery_extra/waveguide")) # looks for modules in the correct directory
using waveguide_debug
using Test

println("--------------------------")
println("--------------------------")
println("")
println("  This tests a reference implementation of the WEP")
println("  You need to download the reference MATLAB implementation to run it")
println("")
println("--------------------------")
println("--------------------------")
sleep(1.5)

@testset "WEP_large" begin

@info "This is WEP-test"

delta = 0.1

matlab_debug_eigval_comp_WEP_FD_and_SPMF(105, 35, delta)

debug_WEP_FD_preconditioner(delta)

matlab_debug_WEP_FD(119, 115, delta)

matlab_debug_full_matrix_WEP_FD_SPMF(21, 17, delta)

debug_sqrt_schur(281)

debug_sqrt_derivative()

debug_Mlincomb_FD_WEP(29, 31, delta)

matlab_debug_Schur_WEP_FD(49, 45, delta)

fft_debug_mateq(431, 427, delta)

debug_Sylvester_SMW_WEP(109, 105, delta, 7)

end


text:
using Espresso: matchex

"""
    map(ex, args)

Returns an instance of the map type corresponding to the given expression.

### Input

- `ex`   -- an expression defining the map, in the form of an anonymous function
- `args` -- additional optional arguments

### Output

A map that best matches the given expression.

### Examples

Let us first create a linear map using this macro:

```jldoctest
julia> @map x -> [1 0; 0 0]*x
LinearMap{Int64,Array{Int64,2}}([1 0; 0 0])
```

We can create an affine system as well:

```jldoctest
julia> @map x -> [1 0; 0 0]*x + [2, 0]
AffineMap{Int64,Array{Int64,2},Array{Int64,1}}([1 0; 0 0], [2, 0])
```

Additional arguments can be passed to `@map` using the function-call form, i.e.
separating the arguments by commas, and using parentheses around the macro call.
For exmple, an identity map of dimension 5 can be defined as:

```jldoctest
julia> @map(x -> x, dim=5)
IdentityMap(5)
```

An identity map can alternatively be created by giving a the size of the identity
matrix as `I(n)`, for example:

```jldoctest
julia> @map x -> I(5)*x
IdentityMap(5)
```
"""
macro map(ex, args)
    quote
        local x = $(ex.args)[1]
        local rhs = $(ex.args)[2].args[2]

        # x -> x, dim=...
        MT = IdentityMap
        local pat = Meta.parse("$x")
        local matched = matchex(pat, rhs)
        matched != nothing && return MT($args)

        throw(ArgumentError("unable to match the given expression to a known map type"))
    end
end

macro map(ex)
    quote
        local x = $(ex.args)[1]
        local rhs = $(ex.args)[2].args[2]

        # x -> I(n)*x
        # (this rule is more specific than x -> Ax so it should come before it)
        MT = IdentityMap
        pat = Meta.parse("I(_n) * $x")
        matched = matchex(pat, rhs)
        matched != nothing && return MT(eval(matched[:_n]))

        # x -> Ax
        MT = LinearMap
        local pat = Meta.parse("_A * $x")
        local matched = matchex(pat, rhs)
        matched != nothing && return MT(eval(matched[:_A]))

        # x -> Ax + b
        MT = AffineMap
        pat = Meta.parse("_A * $x + _b")
        matched = matchex(pat, rhs)
        matched != nothing && return MT(eval(matched[:_A]), eval(matched[:_b]))

        throw(ArgumentError("unable to match the given expression to a known map type"))
    end
end


text:
"""
Julia UnetSocket API.
"""
module UnetSockets

using Reexport
@reexport using Fjage
using Dates, Distributed

export UnetSocket, Protocol, Services, Address, Topics, ReservationStatus
export isclosed, getgateway, host, unbind, isbound, connect, disconnect, isconnected, getlocaladdress
export getlocalprotocol, getremoteaddress, getremoteprotocol, settimeout, gettimeout, cancel

include("Messages.jl")

"Well-known addresses."
module Address
  const BROADCAST = 0
end

"Well-known protocol number assignments."
module Protocol
  const DATA = 0
  const RANGING = 1
  const LINK = 2
  const REMOTE = 3
  const MAC = 4
  const ROUTING = 5
  const TRANSPORT = 6
  const ROUTE_MAINTENANCE = 7
  const LINK2 = 8
  const USER = 32
  const MAX = 63
end

"Well-known topics."
module Topics
  const PARAMCHANGE = "org.arl.unet.Topics.PARAMCHANGE"
  const LIFECYCLE = "org.arl.unet.Topics.LIFECYCLE"
end

"Well-known service names."
module Services
  const NODE_INFO = "org.arl.unet.Services.NODE_INFO"
  const ADDRESS_RESOLUTION = "org.arl.unet.Services.ADDRESS_RESOLUTION"
  const DATAGRAM = "org.arl.unet.Services.DATAGRAM"
  const PHYSICAL = "org.arl.unet.Services.PHYSICAL"
  const RANGING = "org.arl.unet.Services.RANGING"
  const BASEBAND = "org.arl.unet.Services.BASEBAND"
  const LINK = "org.arl.unet.Services.LINK"
  const MAC = "org.arl.unet.Services.MAC"
  const ROUTING = "org.arl.unet.Services.ROUTING"
  const ROUTE_MAINTENANCE = "org.arl.unet.Services.ROUTE_MAINTENANCE"
  const TRANSPORT = "org.arl.unet.Services.TRANSPORT"
  const REMOTE = "org.arl.unet.Services.REMOTE"
  const STATE_MANAGER = "org.arl.unet.Services.STATE_MANAGER"
  const SCHEDULE = "org.arl.unet.Services.SCHEDULE"
  const SHELL = "org.arl.fjage.shell.Services.SHELL"
end

"Status indicator for a particular request during the reservation process."
module ReservationStatus
  const START = 0             # Start of channel reservation for a reservation request
  const END = 1               # End of channel reservation for a reservation request
  const FAILURE = 2           # Failure to reserve channel for a reservation request
  const CANCEL = 3            # Cancel channel reservation for a reservation request
  const REQUEST = 4           # Request information from a client agent for a reservation request
end

"""
    sock = UnetSocket(host, port)

Open a new UnetSocket via TCP/IP to communicate with UnetStack.
"""
mutable struct UnetSocket
  gw::Union{Gateway, Nothing}
  localprotocol::Integer
  remoteaddress::Integer
  remoteprotocol::Integer
  timeout::Integer
  provider::Union{AgentID, Nothing}
  waiting::Bool
  function UnetSocket(host::String, port::Integer)
    sock = new(Gateway(host, port), -1, -1, 0, -1, nothing, false)
    @async begin
      alist = agentsforservice(sock.gw, Services.DATAGRAM)
      for a in alist
        subscribe(sock.gw, a)
      end
    end
    return sock
  end
end

Base.show(io::IO, sock::UnetSocket) = Base.show(io, sock.gw)

"Close the socket."
function close(sock::UnetSocket)
  Base.close(sock.gw)
  sock.gw = nothing
end

"Checks if a socket is closed."
function isclosed(sock::UnetSocket)
  return sock.gw == nothing
end

"Binds a socket to listen to a specific protocol datagrams."
function bind(sock::UnetSocket, protocol::Integer)
  if protocol == Protocol.DATA || (protocol >= Protocol.USER && protocol <= Protocol.MAX)
    sock.localprotocol = protocol
    return true
  end
  return false
end

"Unbinds a socket so that it listens to all unreserved protocols."
function unbind(sock::UnetSocket)
  sock.localprotocol = -1
end

"Checks if a socket is bound."
function isbound(sock::UnetSocket)
  return sock.localprotocol >= 0
end

"""
Sets the default destination address and destination protocol number for
datagrams sent using this socket. The defaults can be overridden for specific
send() calls.
"""
function connect(sock::UnetSocket, to::Integer, protocol::Integer)
  if to >= 0 && (protocol == Protocol.DATA || (protocol >= Protocol.USER && protocol <= Protocol.MAX))
    sock.remoteaddress = to
    sock.remoteprotocol = protocol
    return true
  end
  return false
end

"Resets the default destination address to undefined, and the default protocol number to Protocol.DATA."
function disconnect(sock::UnetSocket)
  sock.remoteaddress = -1
  sock.remoteprotocol = 0
end

"Checks if a socket is connected, i.e., has a default destination address and protocol number."
function isconnected(sock::UnetSocket)
  return sock.remoteaddress >= 0
end

"Gets the local node address."
function getlocaladdress(sock::UnetSocket)
  if sock.gw == nothing
    return -1
  end
  nodeinfo = agentforservice(sock.gw, Services.NODE_INFO)
  if nodeinfo == nothing
    return -1
  end
  req = ParameterReq(recipient=nodeinfo, param="address")
  rsp = request(sock.gw, req)
  if rsp == nothing
    return -1
  end
  return rsp.value
end

"Gets the protocol number that the socket is bound to."
function getlocalprotocol(sock::UnetSocket)
  return sock.localprotocol
end

"Gets the default destination node address for a connected socket."
function getremoteaddress(sock::UnetSocket)
  return sock.remoteaddress
end

"Gets the default transmission protocol number."
function getremoteprotocol(sock::UnetSocket)
  return sock.remoteprotocol
end

"""
Sets the timeout for datagram reception. The default timeout is infinite.
i.e., the receive() call blocks forever. A timeout of 0 means the
receive() call is non-blocking. Negative timeout is considered infinite.
"""
function settimeout(sock::UnetSocket, ms::Integer)
  if ms < 0
    ms = -1
  end
  sock.timeout = ms
end

"Gets the timeout for datagram reception."
function gettimeout(sock::UnetSocket)
  return sock.timeout
end

"Transmits a datagram to the default node address using the default protocol."
function Fjage.send(sock::UnetSocket, data)
  if sock.remoteaddress < 0
    return false
  end
  return send(sock, DatagramReq(to=sock.remoteaddress, protocol=sock.remoteprotocol, data=Vector{UInt8}(data)))
end

"Transmits a datagram to the specified node address using the default protocol."
function Fjage.send(sock::UnetSocket, data, to::Integer)
  return send(sock, DatagramReq(to=to, protocol=sock.remoteprotocol, data=Vector{UInt8}(data)))
end

"Transmits a datagram to the specified node address using the specified protocol."
function Fjage.send(sock::UnetSocket, data, to::Integer, protocol::Integer)
  return send(sock, DatagramReq(to=to, protocol=protocol, data=Vector{UInt8}(data)))
end

"Transmits a datagram to the specified node address using the specified protocol."
function Fjage.send(sock::UnetSocket, req::Message)
  if sock.gw == nothing
    return false
  end
  if !isa(req, DatagramReq)
    return false
  end
  protocol = req.protocol
  protocol == nothing && (protocol = Protocol.USER)
  if protocol != Protocol.DATA && (protocol < Protocol.USER || protocol > Protocol.MAX)
    return false
  end
  if req.recipient == nothing
    (sock.provider == nothing) && (sock.provider = agentforservice(sock.gw, Services.TRANSPORT))
    (sock.provider == nothing) && (sock.provider = agentforservice(sock.gw, Services.ROUTING))
    (sock.provider == nothing) && (sock.provider = agentforservice(sock.gw, Services.LINK))
    (sock.provider == nothing) && (sock.provider = agentforservice(sock.gw, Services.PHYSICAL))
    (sock.provider == nothing) && (sock.provider = agentforservice(sock.gw, Services.DATAGRAM))
    (sock.provider == nothing) && (return false)
    req.recipient = sock.provider
  end
  rsp = request(sock.gw, req)
  return rsp != nothing && rsp.performative == Performative.AGREE
end

"""
Receives a datagram sent to the local node and the bound protocol number. If the
socket is unbound, then datagrams with all unreserved protocols are received. Any
broadcast datagrams are also received.

This call blocks until a datagram is available, the socket timeout is reached,
or until cancel() is called.
"""
function Fjage.receive(sock::UnetSocket)
  if sock.gw == nothing
    return nothing
  end
  t0 = now()
  while sock.timeout <= 0 || (now()-t0).value < sock.timeout
    sock.waiting = true
    ntf = receive(sock.gw, DatagramNtf, sock.timeout)
    sock.waiting = false
    if ntf == nothing
      return nothing
    end
    if isa(ntf, DatagramNtf)
      p = ntf.protocol
      if p == Protocol.DATA || p >= Protocol.USER
        if sock.localprotocol < 0 || sock.localprotocol == p
          return ntf
        end
      end
    end
  end
  return nothing
end

"Cancels an ongoing blocking receive()."
function cancel(sock::UnetSocket)
  if sock.waiting
    push!(sock.gw.queue, nothing)
  end
end

"Gets a Gateway to provide low-level access to UnetStack."
function getgateway(sock::UnetSocket)
  return sock.gw
end

"Gets an AgentID providing a specified service for low-level access to UnetStack."
function Fjage.agentforservice(sock::UnetSocket, svc::String)
  if sock.gw == nothing
    return nothing
  else
    return agentforservice(sock.gw, svc)
  end
end

"Gets a list of AgentIDs providing a specified service for low-level access to UnetStack."
function Fjage.agentsforservice(sock::UnetSocket, svc::String)
  if sock.gw == nothing
    return []
  else
    return agentsforservice(sock.gw, svc)
  end
end

"Gets a named AgentID for low-level access to UnetStack."
function Fjage.agent(sock::UnetSocket, name::String)
  if sock.gw == nothing
    return nothing
  else
    return agent(sock.gw, name)
  end
end

"""
    address = host(sock, name)

Resolve node name to node address.
"""
function host(sock::UnetSocket, name::String)
  if sock.gw == nothing
    return nothing
  end
  arp = agentforservice(sock.gw, Services.ADDRESS_RESOLUTION)
  if arp == nothing
    return nothing
  end
  rsp = arp << AddressResolutionReq(name=name)
  if rsp == nothing
    return nothing
  end
  return rsp.address
end

# Base functions to add local methods
Base.close(sock::UnetSocket) = close(sock)
Base.bind(sock::UnetSocket, protocol::Integer) = bind(sock, protocol)

end


text:
module LazDDA
using HDF5,PyPlot,Statistics,LsqFit,PyCall,FFTW,StatsBase,Images
using LazCore,LazType,LazThermal,LazIO,LazThermal_Kritsuk
using LazRHT_investigation,LazCFA

	##############################################################################
	#
	# Copyright (c) 2020
	# Ka Ho Yuen, Alex Lazarian and Dmitri Pogosyan
	# All Rights Reserved.
	#
	# This program is free software: you can redistribute it and/or modify
	# it under the terms of the GNU General Public License as published by
	# the Free Software Foundation, either version 3 of the License, or
	# (at your option) any later version.
	#
	# This program is distributed in the hope that it will be useful,
	# but WITHOUT ANY WARRANTY; without even the implied warranty of
	# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	# GNU General Public License for more details.
	# You should have received a copy of the GNU General Public License
	# along with this program.  If not, see <https://www.gnu.org/licenses/>.
	#
	##############################################################################

	#==
	 KH: This is v1.3+ only
	==#

	# Nickname "Dmitri Decomposition Algorithm", referring to the algorithm me and Dmitri came with in April 2020
	# when trying to decompose the velocity channel map into densities and velocities contribution.

	export pdec,DDA,DDA_dispersion,vbin

	function re(p::Cube;dims=3,s=size(p)[1:2])
		I=sum(p,dims=dims);
		I=reshape(I,s);
		return I
	end

	function channelproduct(A::Mat,B::Mat)
	    return mean((A.-mean(A)).*(B.-mean(B)))
	end

	function vbin(v)
		vv=zeros(2*length(v)-1);
		vv[1:2:end]=v;
		vv[2:2:end]=0.5.*(v[1:end-1].+v[2:end]);
		return vv
	end

	function ncp(A::Mat,B::Mat)
	    return channelproduct(A,B)./std(A)./std(B)
	end

	function normalized_channel(A::Mat)
	    return (A.-mean(A))./std(A)
	end

	function pd(p::Mat,I::Mat)
	    Ix=normalized_channel(I);
	    return channelproduct(p,Ix).*Ix
	end

	function pv(p::Mat,I::Mat)
	    return p.-pd(p,I)
	end

	function pdec(p::Mat,I::Mat)
	    return pd(p,I),pv(p,I)
	end

	function DDA(p::Cube)
		nx,ny,nv=size(p);
		pdd=zeros(size(p));
		pvv=zeros(size(p));
		I=re(p)
		for k in 1:nv
			pdd[:,:,k]=pd(p[:,:,k],I)
			pvv[:,:,k]=pv(p[:,:,k],I)
		end
		return pdd,pvv
	end

	function DDA_dispersion(p::Cube,vv)
		pd,pv=DDA(p)
		aaaa=zeros(0)
		bbbb=zeros(0)
		cccc=zeros(0)
		for i in 1:length(vv)
		    push!(aaaa,vv[i])
		    push!(bbbb,std(pd[:,:,i]))
		    push!(cccc,std(pv[:,:,i]))
		end
		return aaaa,bbbb,cccc
	end

end # LazDDA


text:
############################################################################################
#import
include("diagnostics.jl")
include("tune.jl")

############################################################################################
#Export


text:
using Roots, Dierckx, NLopt, CoolProp, PyPlot, LsqFit
import JSON
include("rel_perms.jl")
include("frac_flow_funcs.jl")

# read the data files
lowsal=JSON.parsefile("../juliaphreeqc/data/Lowsal_data.json")
paper_title = "fathi2010"
exp_type = "coreFlood"
exp_name = ["SF20-a", "SF20-b"] # formation water and sea water
recovery_data = "coreFloodData" # -> time(day), recovery
core_flood_condition = lowsal["fathi2010"]["coreFlood"]["SF20-a"]
L_core = core_flood_condition["L(mm)"]/1000
D_core = core_flood_condition["D(mm)"]/1000
poros  = core_flood_condition["porosity"]/100
swi     = core_flood_condition["swi"]/100
soi     = core_flood_condition["soi"]/100 # 1-swi
T_c     = core_flood_condition["T(C)"]
oil_type = core_flood_condition["oil"]
mu_oil_25 = lowsal["fathi2010"]["oil"][oil_type]["viscosity(cp)"]/100 # Pa.s
t_day = lowsal["fathi2010"]["coreFloodData"]["SF20-a"]["time(day)"]
R_oil = lowsal["fathi2010"]["coreFloodData"]["SF20-a"]["recovery"]/100
injection_rate_pv_day = core_flood_condition["injectionRate(PV/day)"]
recovery_final = core_flood_condition["recovery(%OOIP)"]/100

so_final = (1-recovery_final)*soi # final oil saturation
pv_inj = t_day[end]*injection_rate_pv_day
t_sec = t_day*24*3600 # s
# conversion of the input data to the right unit
A_core  = π*D_core^2/4
PV_core = L_core*A_core*poros
u_inj   = injection_rate_pv_day/(24*3600)*PV_core/A_core # Darcy velocity
T_K = T_c + 273.15
perm_ave = 0.002e-12 # m^2 permeability
mu_water = PropsSI("viscosity", "T", T_K, "P", 5e5, "H2O")
# My formulation for oil viscosity at higher Temperature
# source: P.Daučík et al., Temperature Dependence of the Viscosity of Hydrocarbon Fractions
# Acta Chimica Slovaca, Vol.1, No. 1, 2008, 43 – 57
# log(mu/mu0) = B(1/T - 1/T0)
# find B for n-Dodecane, and use it for extrapolation
mu_1   = PropsSI("viscosity", "T", 300, "P", 5e5, "n-Dodecane")
mu_2   = PropsSI("viscosity", "T", T_K, "P", 5e5, "n-Dodecane")
B_visc = log(mu_1/mu_2)/(1/300-1/T_K)
mu_oil = exp(B_visc*(1/T_K-1/(25+273.15)))*mu_oil_25 # Pa.s

# Corrections to make the data consistent;
# The recovery line must be a straight line at the begining of the experiment.
# I shift the time axis a bit to make it a straight line:
# plot(t_sec, R_oil, "-o")
(a,b) = linreg(t_sec[2:6], R_oil[2:6])
# 0 = a + b*t_shift
t_shift = -a/b
t_sec_cor = [0; t_sec[2:end]-t_shift]
# plot(t_sec_cor, R_oil, "-o")

# define the objective functions
# (xt_shock_ww, sw_shock_ww, xt_prf_ww, sw_prf_ww, t_ww, p_inj_ww, R_oil_ww)=
# frac_flow_wf(muw=mu_water, muo=mu_oil, ut=inj_vel, phi=poros_res,
#   k=perm_res, swc=swc_ww, sor=sor_ww, kro0=kro0_ww, no=no_ww, krw0=krw0_ww,
#   nw=nw_ww, sw0=swc_ww, sw_inj=1.0, L=L_res, pv_inj=inj_pv)

# define the objective functions
# param = [sor, swc, kro0, krw0, no, nw]
function rec_fact(param, muw, muo, ut, phi, k, sw0, sw_inj, L, pv_inj, t_exp, R_exp)
  (sor, swc, kro0, krw0, no, nw) = param
  (xt_shock, sw_shock, xt_prf, sw_prf, t, p_inj, R)=
  frac_flow_wf(muw=muw, muo=muo, ut=ut, phi=phi,
    k=k, swc=swc, sor=sor, kro0=kro0, no=no, krw0=krw0,
    nw=nw, sw0=sw0, sw_inj=sw_inj, L=L, pv_inj=pv_inj)
  R_int = Spline1D(t, R, k=1)
  return R_int(t_exp)
end

model = (t_exp, param) -> rec_fact(param, mu_water, mu_oil, u_inj, poros, perm_ave, swi, 1.0, L_core, pv_inj, t_sec_cor, R_oil)

function error_calc(param, muw, muo, ut, phi, k, sw0, sw_inj, L, pv_inj, t_exp, R_exp)
  (sor, swc, kro0, krw0, no, nw) = param
  (xt_shock, sw_shock, xt_prf, sw_prf, t, p_inj, R)=
  frac_flow_wf(muw=muw, muo=muo, ut=ut, phi=phi,
    k=k, swc=swc, sor=sor, kro0=kro0, no=no, krw0=krw0,
    nw=nw, sw0=sw0, sw_inj=sw_inj, L=L, pv_inj=pv_inj)
  R_int = Spline1D(t, R, k=1)
  error_vals = R_int(t_exp) - R_exp
  return error_vals
end

function error_calculation(param, param_ind, muw, muo, ut, phi, k, sw0, sw_inj, L, pv_inj,
  t_exp, R_exp; grad_calc=true)
  eps1 = 1e-8
  error_vals = error_calc(param, muw, muo, ut, phi, k, sw0, sw_inj, L, pv_inj, t_exp, R_exp)
  if grad_calc
    jac  = zeros(length(error_vals), length(param_ind))
  else
    jac = zeros(0)
  end
  if grad_calc
    for j in eachindex(param_ind)
      param2 = copy(param)
      param2[param_ind[j]]+=eps1
      error_val2 = error_calc(param2, muw, muo, ut, phi, k, sw0, sw_inj, L, pv_inj, t_exp, R_exp)
      jac[:, j] = (error_val2 - error_vals) / eps1
    end
  end
  return error_vals, jac
end

function objective_function(param::Vector{Float64}, grad::Vector{Float64}, param_all,
  param_ind, muw, muo, ut, phi, k, sw0, sw_inj, L, pv_inj, t_exp, R_exp; w=1.0)
  param2 = copy(param_all)
  param2[param_ind] = param
  grad_calc = false
  if length(grad)>0
    grad_calc = true
  end
  error_val, jac = error_calculation(param2, param_ind, muw, muo, ut, phi, k, sw0,
    sw_inj, L, pv_inj, t_exp, R_exp, grad_calc = grad_calc)
  obj_func_value = sum(abs2, w.*error_val)
  if length(grad)>0
      grad[:] = 2.0*sum(w.*error_val.*jac, 1)
  end
  return obj_func_value
end

function plot_results(param, muw, muo, ut, phi, k, sw0, sw_inj, L, pv_inj, t_exp, R_exp)
  (sor, swc, kro0, krw0, no, nw) = param
  (xt_shock, sw_shock, xt_prf, sw_prf, t, p_inj, R)=
  frac_flow_wf(muw=muw, muo=muo, ut=ut, phi=phi,
    k=k, swc=swc, sor=sor, kro0=kro0, no=no, krw0=krw0,
    nw=nw, sw0=sw0, sw_inj=sw_inj, L=L, pv_inj=pv_inj)
  R_int = Spline1D(t, R, k=1)
  error_vals = R_int(t_exp) - R_exp
  plot(t_exp, R_int(t_exp), t_exp, R_exp, "o")
end

plot_quick = param -> plot_results(param, mu_water, mu_oil, u_inj, poros, perm_ave, swi, 1.0, L_core, pv_inj, t_sec_cor, R_oil)
# Optimization problem:

# parameter values
# [sor, swc, kro0, krw0, no, nw]
param_all = [so_final, 0.25, 0.7, 0.3, 2, 2.4]
param_ind = [1, 2,3,4,5,6]
w = ones(length(R_oil))
w[end]=5
w[end-1]=15
w[7] = 3
obj_fun = (param, grad) -> objective_function(param::Vector{Float64}, grad::Vector{Float64}, param_all,
  param_ind, mu_water, mu_oil, u_inj, poros, perm_ave, swi, 1.0, L_core, pv_inj, t_sec_cor, R_oil, w=w)

# x_init = [0.2, 0.8, 0.4, 2, 2]
# x_lb = [0.1, 0.1, 0.1, 0.5, 0.5]
# x_ub = [0.5, 1.0, 1.0, 5.0, 5.0]

x_init = copy(param_all)
x_lb = [0.1, 0.1, 0.05, 0.05, 1.0, 1.0]
x_ub = [0.5, 0.5, 1.0, 1.0, 4.0, 4.0]

fit = curve_fit(model, t_sec_cor, R_oil, w, x_init, lower = x_lb, upper = x_ub)
#
# # algorithms
# # :LD_MMA
# # :LN_COBYLA
# # :LD_LBFGS
# # :GN_DIRECT
# # :GN_DIRECT_L
# opt_alg=:LN_COBYLA
#
# opt1 = Opt(opt_alg, length(x_init)) # choose the algorithm
# lower_bounds!(opt1, x_lb)
# upper_bounds!(opt1, x_ub)
# ftol_rel!(opt1, 1e-10)
# ftol_abs!(opt1, 1e-10)
#
# min_objective!(opt1, obj_fun)
# (fObjOpt, paramOpt, flag) = optimize(opt1, x_init)
#
# param_plot = copy(param_all)
# param_plot[param_ind] = paramOpt
plot_quick(fit.param)

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng


text:
# This file is an example for the chemfiles library
# Any copyright is dedicated to the Public Domain.
# http://creativecommons.org/publicdomain/zero/1.0/
#!/usr/bin/env julia
using Chemfiles

trajectory = Trajectory("input.arc")
output = Trajectory("output.pdb", 'w')

selection = Selection("name Zn or name N")

for frame in trajectory
    to_remove = evaluate(selection, frame)
    for i in reverse(sort(to_remove))
        remove_atom!(frame, i)
    end
    write(output, frame)
end

close(trajectory)
# When running on the REPL, closing the trajectory
# is required to flush buffered output and finish writing
close(output)


text:
using Distributions
using PyPlot
using LsqFit
using CSV, DataFrames, DataFramesMeta
using StatsBase

""" Functions for results of Fig. 2E """
# Run:
# mean_thrs_delay_1, mean_thrs_delay_3 = SinergiafMRI_datafit.get_state_visits_bootstrapped()
function get_state_visits_bootstrapped(;
                                        # filestr = "_Sim_",
                                        filestr = ""
                                        )
    # Real:
    data_path = "./mcmc_rl_fit/projects/fmri/data/SARSPEICZVG_all_fMRI.csv"
    # Simulated:
    # data_path = "./mcmc_rl_fit/projects/fmri/data_sim/SARSPEICZVG_all_fMRI_Sim_exporder_SurpAC_run2.csv"

    max_episodes_per_subj = 40
    (all_data, all_stats) = RL_Fit.load_SARSPEI(data_path;max_episodes_per_subj=max_episodes_per_subj)
    all_subject_ids = all_stats.all_subject_ids
    println("all_subject_ids: $all_subject_ids")

    nrOfBootstrapRounds = 200

    threshold_delay_matrices_1 = Array{Array{Float64,1}, 1}(undef, nrOfBootstrapRounds)
    [threshold_delay_matrices_1[i] = zeros(Float64, 7) for i in 1:nrOfBootstrapRounds]
    threshold_delay_matrices_3 = Array{Array{Float64,1}, 1}(undef, nrOfBootstrapRounds)
    [threshold_delay_matrices_3[i] = zeros(Float64, 7) for i in 1:nrOfBootstrapRounds]
    #
    # threshold_delay_matrices_nofitting_1 = Array{Array{Float64,1}, 1}(undef, nrOfBootstrapRounds)
    # [threshold_delay_matrices_nofitting_1[i] = zeros(Float64, 7) for i in 1:nrOfBootstrapRounds]
    # threshold_delay_matrices_nofitting_3 = Array{Array{Float64,1}, 1}(undef, nrOfBootstrapRounds)
    # [threshold_delay_matrices_nofitting_3[i] = zeros(Float64, 7) for i in 1:nrOfBootstrapRounds]

    graph_1_subjIDs = [4, 5, 6, 8, 10, 12, 15, 16, 19, 21]
    graph_3_subjIDs = [1, 2, 3, 7, 9, 11, 13, 14, 17, 18, 20]
    for i in 1:nrOfBootstrapRounds
        # Sample with replacement

        bootstrapped_subject_ids_graph_1 = sample(graph_1_subjIDs, length(graph_1_subjIDs), replace=true)
        bootstrapped_subject_ids_graph_3 = sample(graph_3_subjIDs, length(graph_3_subjIDs), replace=true)
        bootstrapped_subject_ids = vcat(bootstrapped_subject_ids_graph_1, bootstrapped_subject_ids_graph_3)

        # Get graph_visits
        graph_1_visits, graph_3_visits = get_state_visits_singleround(all_data,
                        bootstrapped_subject_ids, max_episodes_per_subj = max_episodes_per_subj)

        threshold_delay_matrices_1[i] = get_state_fittedthreshold(1, graph_1_visits)
        threshold_delay_matrices_3[i] = get_state_fittedthreshold(3, graph_3_visits)

        # threshold_delay_matrices_nofitting_1[i] = get_state_threshold_nofitting(1, graph_1_visits)
        # threshold_delay_matrices_nofitting_3[i] = get_state_threshold_nofitting(3, graph_3_visits)

    end

    mean_thrs_delay_1, std_thrs_delay_1, stderror_thrs_delay_1 =  get_bootstrapped_stats(threshold_delay_matrices_1)
    mean_thrs_delay_3, std_thrs_delay_3, stderror_thrs_delay_3 =  get_bootstrapped_stats(threshold_delay_matrices_3)

    write_bootstrapped_results_csv(1, mean_thrs_delay_1, std_thrs_delay_1, stderror_thrs_delay_1, filestr = filestr)
    write_bootstrapped_results_csv(3, mean_thrs_delay_3, std_thrs_delay_3, stderror_thrs_delay_3, filestr = filestr)

    # ------------------------------------------
    # mean_thrs_delay_nofitting_1, std_thrs_delay_nofitting_1, stderror_thrs_delay_nofitting_1 =  get_bootstrapped_stats(threshold_delay_matrices_nofitting_1)
    # mean_thrs_delay_nofitting_3, std_thrs_delay_nofitting_3, stderror_thrs_delay_nofitting_3 =  get_bootstrapped_stats(threshold_delay_matrices_nofitting_3)
    #
    # write_bootstrapped_results_csv(1, mean_thrs_delay_nofitting_1, std_thrs_delay_nofitting_1, stderror_thrs_delay_nofitting_1, filestr=filestr*"nofitting_")
    # write_bootstrapped_results_csv(3, mean_thrs_delay_nofitting_3, std_thrs_delay_nofitting_3, stderror_thrs_delay_nofitting_3, filestr=filestr*"nofitting_")

    mean_thrs_delay_1, mean_thrs_delay_3#,  mean_thrs_delay_nofitting_1, mean_thrs_delay_nofitting_3
end

""" Similar to get_state_visits() """
function get_state_visits_singleround(all_data, subject_ids; max_episodes_per_subj = 40)

    graph_1_visits = zeros(Int64, max_episodes_per_subj,7, 2)
    graph_3_visits = zeros(Int64, max_episodes_per_subj,7, 2)
    # @show subject_ids
    for subj_id in subject_ids
        # @show subj_id
        subj_data = RL_Fit.filter_by_person(all_data, subj_id)
        graph_id = subj_data[1,Int(RL_Fit.c_graph_version)]
        all_episodes = unique(subj_data[:,Int(RL_Fit.c_episode)])
        trajectory_count_per_state = zeros(Int64, 7) #count for each state individually how many times it participated in a rewarded episode
        # per subject, do not count (sum), just set to 0 or 1.
        state_action_indicator = zeros(Int64, max_episodes_per_subj, 7, 2)

        for episode in all_episodes
            for s in 2:7
                # check if state s has "contributed" to the current episode. if yes, count this trajectory.
                if any( (subj_data[:,Int(RL_Fit.c_episode)] .== episode) .& (subj_data[:,Int(RL_Fit.c_state)] .== s) )
                    trajectory_count_per_state[s] += 1
                    # now count the value of actions. distinguish two cases:
                    # for some states, the two actions are equally valid. In this case count the action
                    # for other states, there is a correct (1) and a wrong(-1) action
                    # action correct:
                    if any( (subj_data[:,Int(RL_Fit.c_episode)] .== episode) .& (subj_data[:,Int(RL_Fit.c_state)] .== s) .& (subj_data[:,Int(RL_Fit.c_action_value)] .== 0))
                        # action A or action B?
                        if any( (subj_data[:,Int(RL_Fit.c_episode)] .== episode) .& (subj_data[:,Int(RL_Fit.c_state)] .== s) .& (subj_data[:,Int(RL_Fit.c_action)] .== 1))
                            # A
                            state_action_indicator[trajectory_count_per_state[s], s, 1] = 1 # do not +=1
                        end
                        # note: it's not an ELSE here. we check for both and count both options at most once
                        if any( (subj_data[:,Int(RL_Fit.c_episode)] .== episode) .& (subj_data[:,Int(RL_Fit.c_state)] .== s) .& (subj_data[:,Int(RL_Fit.c_action)] .== 2))
                            # B
                            state_action_indicator[trajectory_count_per_state[s], s, 2] = 1 # do not +=1
                        end
                    else
                        # action correct or wrong?
                        if any( (subj_data[:,Int(RL_Fit.c_episode)] .== episode) .& (subj_data[:,Int(RL_Fit.c_state)] .== s) .& (subj_data[:,Int(RL_Fit.c_action_value)] .== +1))
                            # correct
                            state_action_indicator[trajectory_count_per_state[s], s, 2] = 1 # do not +=1
                        end
                        # note: it's not an ELSE here. we check for both and count both options at most once
                        if any( (subj_data[:,Int(RL_Fit.c_episode)] .== episode) .& (subj_data[:,Int(RL_Fit.c_state)] .== s) .& (subj_data[:,Int(RL_Fit.c_action_value)] .== -1))
                            # wrong
                            state_action_indicator[trajectory_count_per_state[s], s, 1] = 1 # do not +=1
                        end
                    end
                end
            end # end all states
        end # end all episodes
        # @show state_action_indicator
        if graph_id == 1
            graph_1_visits += state_action_indicator
        elseif graph_id == 3
            graph_3_visits += state_action_indicator
        else
            error("unknown graph id: $graph_id")
        end
    end # end all subject
    graph_1_visits, graph_3_visits
end

function get_state_fittedthreshold(graph_id::Int,
                            graph_state_counts;
                            learning_threshold=0.8,
                            min_visits_requirement = 2,
                            savepath = "/Users/vasia/Desktop/")
    dfresults = []
    graph_info = get_graph_info()
    states, minimal_dist, state_txt, graph = get_this_graph_info(graph_info, graph_id)

    init_vals_exp = [0.5, 0.]
    threshold_delay_matrix = zeros(length(states)+1)
    for s = states
        # println(s)
        counts_per_episode_actionwrong = vec(graph_state_counts[:,s,1])
        counts_per_episode_actioncorrect = vec(graph_state_counts[:,s,2])
        summed_counts_per_episode = counts_per_episode_actionwrong.+counts_per_episode_actioncorrect
        episodes_min_visits = findall(summed_counts_per_episode .>= min_visits_requirement)
        # @show episodes_min_visits
        ratio = counts_per_episode_actioncorrect[episodes_min_visits] ./ summed_counts_per_episode[episodes_min_visits]
        weight = (summed_counts_per_episode[episodes_min_visits])
        nr_x = length(ratio)

        x_data = 1:nr_x
        y_data = ratio
        # markersize = 1+2*sqrt(weight[i]),
        threshold_delay = 0.
        weight = Float64.(weight)
        if minimal_dist[1,s] != minimal_dist[2,s]
            fit = curve_fit(model_exp, x_data, y_data, weight, init_vals_exp) # Weighted cost function
            y_hat = model_exp(x_data, fit.param)
            threshold_delay = inv_model_exp(learning_threshold, fit.param)
        end
        if threshold_delay < 0.
            if fit.param[2] > 0.
                threshold_delay = 0.
            else
                threshold_delay = Float64(episodes_min_visits[end])
            end
            # @show threshold_delay
            # println("----------")
        end
        if threshold_delay > Float64(episodes_min_visits[end])
            threshold_delay = Float64(episodes_min_visits[end])
        end
        threshold_delay_matrix[s] = threshold_delay
    end
    threshold_delay_matrix
end

# NOTE: NOT USED!!!
function get_state_threshold_nofitting(graph_id::Int,
                            graph_state_counts;
                            learning_threshold=0.8,
                            min_visits_requirement = 2,
                            savepath = "/Users/vasia/Desktop/")
    dfresults = []
    graph_info = get_graph_info()
    states, minimal_dist, state_txt, graph = get_this_graph_info(graph_info, graph_id)

    init_vals_exp = [0.5, 0.]
    threshold_delay_matrix = zeros(length(states)+1)
    for s = states
        # println(s)
        counts_per_episode_actionwrong = vec(graph_state_counts[:,s,1])
        counts_per_episode_actioncorrect = vec(graph_state_counts[:,s,2])
        summed_counts_per_episode = counts_per_episode_actionwrong.+counts_per_episode_actioncorrect
        episodes_min_visits = findall(summed_counts_per_episode .>= min_visits_requirement)
        # @show episodes_min_visits
        ratio = counts_per_episode_actioncorrect[episodes_min_visits] ./ summed_counts_per_episode[episodes_min_visits]
        weight = (summed_counts_per_episode[episodes_min_visits])
        nr_x = length(ratio)

        # markersize = 1+2*sqrt(weight[i]),
        threshold_delay = 0.
        weight = Float64.(weight)
        if minimal_dist[1,s] != minimal_dist[2,s]
            threshold_delay = findfirst(x-> x.>=0.8, ratio)
        end
        if isnothing(threshold_delay)
            threshold_delay = Float64(episodes_min_visits[end])
        end
        threshold_delay_matrix[s] = threshold_delay
    end
    threshold_delay_matrix
end

function get_bootstrapped_stats(threshold_delay_matrices)
    thrs_delay_cat = hcat(threshold_delay_matrices...)
    mean_thrs_delay = mean(thrs_delay_cat, dims=2)
    std_thrs_delay = std(thrs_delay_cat, dims=2)
    @show size(threshold_delay_matrices,1)
    stderror_thrs_delay = std(thrs_delay_cat, dims=2) ./ sqrt(size(threshold_delay_matrices,1))
    mean_thrs_delay, std_thrs_delay, stderror_thrs_delay
end

function write_bootstrapped_results_csv(graph_id::Int64,
                                        mean_thrs_delay,
                                        std_thrs_delay,
                                        stderror_thrs_delay;
                                        filestr="",
                                        # filestr = "_Sim_",
                                        # savepath = "./"
                                        )
    savepath = makehomesavepath("learningcurvebars_bootstrap")
    mkdir(savepath)

    graph_info = get_graph_info()
    states, minimal_dist, state_txt, graph = get_this_graph_info(graph_info, graph_id)

    dfresults = DataFrame()
    for s = states
        min_dist = sort(minimal_dist[:,s] .- 1.)
        dfresults[!, Symbol(string(s)*"_"*state_txt[s]*"_thresholddelay_mean")] = [mean_thrs_delay[s]]
        dfresults[!, Symbol(string(s)*"_"*state_txt[s]*"_thresholddelay_std")] = [std_thrs_delay[s]]
        dfresults[!, Symbol(string(s)*"_"*state_txt[s]*"_thresholddelay_stderror")] = [stderror_thrs_delay[s]]
        dfresults[!, Symbol(string(s)*"_"*state_txt[s]*"_dist1vsdist2")] = [string(Int(min_dist[1]))*"-or-"*string(Int(min_dist[2]))]
        dfresults[!, Symbol(string(s)*"_"*state_txt[s]*"_graphid")] = [graph_id]
    end
    CSV.write(joinpath(savepath, "learningcurvebar_bootstrap_"*filestr*"G" * string(graph_id) *".csv"), dfresults,
        delim = " ");
    # dfresults
end


text:
using Gadfly, Colors, Test

set_default_plot_size(6inch, 3inch)

# This should throw a MethodError, not a stackoverflow
@test_throws MethodError Scale.color_discrete_manual([RGBA(0,0,0,0.1)])

plot([sin, cos], 0, 25, Scale.color_discrete_manual("red", "purple"))


text:
function parse_type(t::Type, supporteddts::Vector{DataType})
	node = ASTNode(:start)
	node.args[:type] = t
	node.args[:supporteddts] = supporteddts
	node
end


text:
@testset "687.longest-univalue-path.jl" begin
    @test longest_univalue_path(TreeNode{Int}([5,4,5,1,1,5])) == 2
    @test longest_univalue_path(TreeNode{Int}([1,4,5,4,4,5])) == 2
end


text:
function ensuresize{T,D}(array::Array{T,D}, dims::Int...)
  @assert(ndims(array) == length(dims), "$(size(array)) is different to $(dims)")

  if size(array) == dims
    return array
  end

  return Array(T, dims...)
end

macro ensuresize(ex)
  array = ex.args[1]
  other = ex.args[2:length(ex.args)]

  expr = Expr(:call, :ensuresize, array, other...)
  return esc(Expr(:(=), array, expr))
end

@doc doc"Ensure that the size of the of the array is at least dims

The inner storage might be preserved
" ->
function ensuresize!{D}(m::DenseArray{D}, dims::UInt...)
end

macro stabilize(ex)
    @assert ex.head == :(=)
    @assert length(ex.args) == 2
    @assert ex.args[1].head == :(::)

    value = ex.args[2]
    valueType = ex.args[1].args[2]

    esc(quote
      if !isa($value, $valueType)
        Base.error("Cannot stabilize type to $valueType for $value")
      end
      $ex
    end)

end

@doc "Fill an array with random number from a uniform distribution" ->
function uniform!{F<:Float}(a::DenseArray{F}, min::F, max::F)
  r::F = max - min
  @inbounds for i = 1:length(a)
    a[i] = rand(F) * r + min
  end
end

@doc "Fill an array with random number from a gaussian distribution" ->
function randn!{F<:Float}(a::DenseArray{F}, mu::F, sigma::F)
  @inbounds for i = 1:length(a)
    a[i] = randn(F) * sigma + mu
  end
end

# ---

function randn!{T,N}(x::DenseArray{T,N})
  @inbounds for i in 1:length(x)
    x[i] = randn()
  end
end


text:
# Activate environment

include("../word_count_helpers.jl")
using FileIO,JSON

# Grab the two arguments that are passed in
# The first is the input file name and the second is the output file name
#inputfile = ARGS[1]
#outputfile = ARGS[2]
inOutFileList = readlines(open(ARGS[1],"r+"))
print(inOutFileList)

# Read in file and clean the text
for line in inOutFileList
    inputfile = split(line)[1]
    outputfile = split(line)[2]

    # Read in file and clean the text
    f = open(inputfile)
    text=cleantext(readlines(f))

    # Count number of times each word appears
    wordCounts = countwords(text)

    # Save the word counts to the passed in output file name
    open(outputfile*".json","w") do f
    	JSON.print(f,wordCounts)
    end
end

max_stars_repo_name: llsc-supercloud/teaching-examples

max_issues_repo_name: llsc-supercloud/teaching-examples

max_forks_repo_name: llsc-supercloud/teaching-examples