text:
function UHF(Alg::UHFa)
    ints = IntegralHelper{Float64}()
    UHF(ints, Alg)
end

function UHF(ints::IntegralHelper{Float64}, Alg::UHFa)
    Fermi.HartreeFock.uhf_header()
    output("Collecting necessary integrals...")
    t = @elapsed begin
        ints["S"]
        ints["T"]
        ints["V"]
        ints["ERI"]
    end
    output("Done in {:10.5f} s", t)

    guess = Options.get("scf_guess")
    if guess == "core"
        Cα, Λ =  RHF_core_guess(ints)
        Cβ = deepcopy(Cα)
    elseif guess == "gwh"
        Cα, Λ = RHF_gwh_guess(ints)
        Cβ = deepcopy(Cα)
    end
    UHF(ints, Cα, Cβ, Λ, Alg)
end

function UHF(ints::IntegralHelper{Float64, <:AbstractERI, AtomicOrbitals}, Cα::AbstractMatrix, Cβ::AbstractMatrix, Λ::AbstractMatrix, Alg::UHFa)
    molecule = ints.molecule
    output(Fermi.string_repr(molecule))

    # Grab options
    maxit = Options.get("scf_max_iter")
    Etol = Options.get("scf_e_conv")
    Dtol = Options.get("scf_max_rms")
    do_diis = Options.get("diis")
    oda = Options.get("oda")
    oda_cutoff = Options.get("oda_cutoff")
    oda_shutoff = Options.get("oda_shutoff")

    # Loop variables
    ite = 1
    E = 0.0
    ΔE = 1.0
    Drms = 1.0
    diis = false
    damp = 0.0
    converged = false

    Nα = molecule.Nα
    Nβ = molecule.Nβ
    S = ints["S"]
    T = ints["T"]
    V = ints["V"]
    ERI = ints["ERI"]
    m = size(S)[1]
    Dα = zeros(Float64, (m,m))
    Dβ = zeros(Float64, (m,m))
    Jα = zeros(Float64, (m,m))
    Jβ = zeros(Float64, (m,m))
    Kα = zeros(Float64, (m,m))
    Kβ = zeros(Float64, (m,m))
    ϵα = zeros(Float64, (m))
    ϵβ = zeros(Float64, (m))
    Fα = zeros(Float64, (m,m))
    Fβ = zeros(Float64, (m,m))
    Dsα = deepcopy(Dα)
    Dsβ = deepcopy(Dβ)
    Fsα = deepcopy(Fα)
    Fsβ = deepcopy(Fβ)

    H = T + V
    Dα = buildD!(Dα, Cα, Nα)
    Dβ = buildD!(Dβ, Cβ, Nβ)
    Dα_old = deepcopy(Dα)
    Dβ_old = deepcopy(Dβ)
    Dsα = deepcopy(Dα)
    Dsβ = deepcopy(Dβ)

    build_fock!(Fα, Fβ, Jα, Jβ, Kα, Kβ, Dα, Dβ, ints)
    output(" Guess Energy {:20.14f}", UHFEnergy(H, Dα, Dβ, Fα, Fβ, molecule.Vnuc))

    output("\n Iter.   {:>15} {:>10} {:>10} {:>8} {:>8} {:>8}", "E[UHF]", "ΔE", "Dᵣₘₛ", "t", "DIIS", "damp")
    output(repeat("-",80))
    if do_diis
        DM = Fermi.DIIS.DIISManager{Float64,Float64}(size=Options.get("ndiis"))
        diis_start = Options.get("diis_start")
    end
    t = @elapsed while ite <= maxit
        t_iter = @elapsed begin
            E_old = E

            # Transform Fock matrices to MO basis
            F̃α = Λ*Fα*Λ
            F̃β = Λ*Fβ*Λ

            # Solve for eigenvalues and eigenvectors
            ϵα, C̃α = LinearAlgebra.eigen(Symmetric(F̃α), sortby=x->x)
            ϵβ, C̃β = LinearAlgebra.eigen(Symmetric(F̃β), sortby=x->x)

            # Transform orbital coefficient matrices to AO basis
            Cα = Λ*C̃α
            Cβ = Λ*C̃β

            # Build density matrices
            buildD!(Dα, Cα, Nα)
            buildD!(Dβ, Cβ, Nβ)

            # Build Fock matrix
            build_fock!(Fα, Fβ, Jα, Jβ, Kα, Kβ, Dα, Dβ, ints)

            # Calculate energy
            E = UHFEnergy(H, Dα, Dβ, Fα, Fβ, molecule.Vnuc)

            # Store vectors for DIIS
            if do_diis
                err_α = transpose(Λ)*(Fα*Dα*S - S*Dα*Fα)*Λ
                err_β = transpose(Λ)*(Fβ*Dβ*S - S*Dβ*Fβ)*Λ
                err_v = vcat(err_α, err_β)
                F_v = vcat(Fα, Fβ)
                push!(DM, F_v, err_v)
            end

            # Branch for ODA vs DIIS convergence aids
            diis = false
            damp = 0.0
            # Use ODA damping?
            if oda && Drms > oda_cutoff && ite < oda_shutoff
                damp = odadamping(Dα, Dsα, Fα, Fsα)
                damp = odadamping(Dβ, Dsβ, Fβ, Fsβ)
            # Or Use DIIS?
            elseif do_diis && ite > diis_start
                diis = true
                F_v = Fermi.DIIS.extrapolate(DM)
                Fα .= F_v[1:m, :]
                Fβ .= F_v[m+1:2m, :]
            end

            # Calculate energy difference, Drms, and check for convergence
            ΔE = E-E_old
            ΔDα = Dα - Dα_old
            ΔDβ = Dβ - Dβ_old
            Drms = (sum(ΔDα.^2)/m^2)^(1/2) + (sum(ΔDβ.^2)/m^2)^(1/2)
            Dα_old .= Dα
            Dβ_old .= Dβ
        end
        output("    {:<3} {:>15.10f} {:>11.3e} {:>11.3e} {:>8.2f} {:>8}    {:5.2f}", ite, E, ΔE, Drms, t_iter, diis, damp)
        ite += 1
        if (abs(ΔE) <= Etol) & (Drms <= Dtol) & (ite > 5)
            converged = true
            break
        end
    end

    nocc = Nα + Nβ  # TODO
    nvir = 2*m - nocc   #    output(repeat("-",80))

    output("    UHF done in {:>5.2f}s", t)
    output("    @Final UHF Energy     {:>20.12f} Eₕ", E)
    output("\n   • Orbitals Summary",)
    output("\n   ⬗ Alpha (α) orbitals")
    output("\n {:>10}   {:>15}   {:>10}", "Orbital", "Energy", "Occupancy")
    for i in eachindex(ϵα)
        output(" {:>10}   {:> 15.10f}   {:>6}", i, ϵα[i], (i ≤ Nα ? "↿" : ""))
    end
    output("\n   ⬗ Beta (β) orbitals")
    output("\n {:>10}   {:>15}   {:>10}", "Orbital", "Energy", "Occupancy")
    for i in eachindex(ϵβ)
        output(" {:>10}   {:> 15.10f}   {:>6}", i, ϵβ[i], (i ≤ Nβ ? "⇂" : ""))
    end
    output("")
    if converged
        output("   ✔  SCF Equations converged 😄")
    else
        output("❗ SCF Equations did not converge in {:>5} iterations ❗", maxit)
    end
    output(repeat("-",80))

    Orbitals = UHFOrbitals(molecule, ints.basis, ϵα, ϵβ, E, Cα, Cβ)
    return UHF(molecule, E, nocc, nvir, Orbitals, ΔE, Drms)
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
using Measurements

5.23 ± 0.14 === 5.23 ± 0.14

(5.23± 0.14) - (5.23 ± 0.14)

(5.23 ± 0.14) / (5.23 ± 0.14)

x = 5.23 ± 0.14
x === x

x - x

x / x

using DifferentialEquations, Measurements, Plots

# Half-life and mean lifetime of radiocarbon, in years
t_12 = 5730 ± 40
τ = t_12 / log(2)

#Setup
u₀ = 1 ± 0
tspan = (0.0, 10000.0)

#Define the problem
radioactivedecay(u,p,t) = - u / τ

#Pass to solver
prob = ODEProblem(radioactivedecay, u₀, tspan)
sol = solve(prob, Tsit5(), reltol = 1e-8)

# Analytic solution
u = exp.(- sol.t / τ)

plot(sol.t, sol.u, label = "Numerical", xlabel = "Years", ylabel = "Fraction of Carbon-14")
plot!(sol.t, u, label = "Analytic")

println("Quantity of carbon-14 after ",  sol.t[11], " years:")
println("Numerical: ", sol[11])
println("Analytic:  ", u[11])

using DifferentialEquations, Measurements, Plots

g = 9.79 ± 0.02; # Gravitational constants
L = 1.00 ± 0.01; # Length of the pendulum

#Initial Conditions
u₀ = [0 ± 0, π / 60 ± 0.01] # Initial speed and initial angle
tspan = (0.0, 6.3)

#Define the problem
function simplependulum(du,u,p,t)
    θ  = u[1]
    dθ = u[2]
    du[1] = dθ
    du[2] = -(g/L)*θ
end

#Pass to solvers
prob = ODEProblem(simplependulum, u₀, tspan)
sol = solve(prob, Tsit5(), reltol = 1e-6)

# Analytic solution
u = u₀[2] .* cos.(sqrt(g / L) .* sol.t)

plot(sol.t, getindex.(sol.u, 2), label = "Numerical")
plot!(sol.t, u, label = "Analytic")

plot(sol.t, getindex.(sol.u, 2) .- u, label = "")

g = 9.79 ± 0.02; # Gravitational constants
L = 1.00 ± 0.01; # Length of the pendulum

#Initial Conditions
u₀ = [0 ± 0, π / 3 ± 0.02] # Initial speed and initial angle
tspan = (0.0, 6.3)

#Define the problem
function simplependulum(du,u,p,t)
    θ  = u[1]
    dθ = u[2]
    du[1] = dθ
    du[2] = -(g/L) * sin(θ)
end

#Pass to solvers
prob = ODEProblem(simplependulum, u₀, tspan)
sol = solve(prob, Tsit5(), reltol = 1e-6)

plot(sol.t, getindex.(sol.u, 2), label = "Numerical")

using DiffEqTutorials
DiffEqTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

ext: jl


text:
module NetCDF
using Compat
using Formatting
using Base.Cartesian
include("netcdf_c.jl")
import Base.show
export NcDim,NcVar,NcFile,ncread,ncread!,ncwrite,nccreate,ncsync,ncinfo,ncclose,ncputatt,NC_BYTE,NC_SHORT,NC_INT,NC_FLOAT,NC_DOUBLE,NC_STRING,ncgetatt,NC_NOWRITE,NC_WRITE,NC_CLOBBER,NC_NOCLOBBER,NC_CLASSIC_MODEL,NC_64BIT_OFFSET,NC_NETCDF4
NC_VERBOSE=false
#Some constants

function __init__()
  global const jltype2nctype=@Compat.Dict(Int8=>NC_BYTE,
                   Int16=>NC_SHORT,
                   Int32=>NC_INT,
                   Int64=>NC_INT64,
                   Float32=>NC_FLOAT,
                   Float64=>NC_DOUBLE,
                   UInt8=>NC_CHAR,
                   ASCIIString=>NC_STRING)

  global const nctype2jltype=@Compat.Dict(NC_BYTE=>Int8,
                    NC_SHORT=>Int16,
                    NC_INT=>Int32,
                    NC_LONG=>Int64,
                    NC_FLOAT=>Float32,
                    NC_DOUBLE=>Float64,
                    NC_CHAR=>UInt8,
                    NC_STRING=>ASCIIString)

  global const nctype2string=@Compat.Dict(NC_BYTE=>"BYTE",
                   NC_SHORT=>"SHORT",
                   NC_INT=>"INT",
                   NC_FLOAT=>"FLOAT",
                   NC_DOUBLE=>"DOUBLE",
                   NC_INT64=>"INT64",
                   NC_STRING=>"STRING")
end

"""The type `NcDim` Represents an NcDim object representing a NetCDF dimension of name `name` with dimension values.
"""
type NcDim
  ncid::Int32
  dimid::Int32
  varid::Int32
  name::ASCIIString
  dimlen::UInt
  vals::AbstractArray
  atts::Dict
end

"""

    NcDim(name::String,dimlength::Integer;values::Union{AbstractArray,Number}=[],atts::Dict{Any,Any}=Dict{Any,Any}())`
This constructor creates an NcDim object with the name `name` and length `dimlength`.
"""
function NcDim(name::AbstractString,dimlength::Integer;values::Union{AbstractArray,Number}=[],atts::Dict{Any,Any}=Dict{Any,Any}())
    (length(values)>0 && length(values)!=dimlength) ? error("Dimension value vector must have the same length as dimlength!") : nothing
    NcDim(-1,-1,-1,ascii(name),dimlength,values,atts)
  end

"""
    NcDim(name::AbstractString,dimlength::Integer;values::Union{AbstractArray,Number}=[],atts::Dict{Any,Any}=Dict{Any,Any}())
This constructor creates an NcDim object with the name `name` and and associated values `values`. Upon creation of the NetCDF file a
dimension variable will be generated and the values be written to this variable. Optionally a Dict of attributes can be supplied.
"""
NcDim(name::AbstractString,values::AbstractArray;atts::Dict=Dict{Any,Any}())=
  NcDim(name,length(values),values=values,atts=atts)
NcDim(name::AbstractString,values::AbstractArray,atts::Dict)=
  NcDim(name,length(values),values=values,atts=atts)

"""
The type `NcVar{T,N}` represents a NetCDF variable. It is a subtype of AbstractArray{T,N}, so normal indexing using `[]`
will work for reding and writing data to and from a NetCDF file. `NcVar` objects can be returned by the `open` function , by
indexing an `NcFIle` object (e.g. `myfile["temperature"]`) or, when creating a new file, constructing it directly.
"""
type NcVar{T,N} <: AbstractArray{T,N}
  ncid::Int32
  varid::Int32
  ndim::Int32
  natts::Int32
  nctype::Int32
  name::ASCIIString
  dimids::Vector{Int32}
  dim::Vector{NcDim}
  atts::Dict
  compress::Int32
end

Base.convert{S,T,N}(::Type{NcVar{T,N}},v::NcVar{S,N})=NcVar{T,N}(v.ncid,v.varid,v.ndim,v.natts,v.nctype,v.name,v.dimids,v.dim,v.atts,v.compress)

"""
    NcVar(name::AbstractString,dimin::Union{NcDim,Array{NcDim,1}};atts::Dict{Any,Any}=Dict{Any,Any}(),t::Union{DataType,Integer}=Float64,compress::Integer=-1)
Here varname is the name of the variable, dimlist an array of type NcDim holding the dimensions associated to the variable, varattributes is a Dict
holding pairs of attribute names and values. t is the data type that should be used for storing the variable. You can either specify a julia type
(Int16, Int32, Float32, Float64) which will be translated to (NC_SHORT, NC_INT, NC_FLOAT, NC_DOUBLE) or directly specify one of the latter list.
You can also set the compression level of the variable by setting compress to a number in the range 1..9 This has only an effect in NetCDF4 files.
"""
#Some additional constructors
function NcVar(name::AbstractString,dimin::Union{NcDim,Array{NcDim,1}};atts::Dict{Any,Any}=Dict{Any,Any}(),t::Union{DataType,Integer}=Float64,compress::Integer=-1)
  dim = isa(dimin,NcDim) ? NcDim[dimin] : dimin
  return NcVar{typeof(t)==DataType ? t : nctype2jltype[t],length(dim)}(-1,-1,length(dim),length(atts), typeof(t)==DataType ? jltype2nctype[t] : t,name,Array(Int32,length(dim)),dim,atts,compress)
end
NcVar(name::AbstractString,dimin::Union{NcDim,Array{NcDim,1}},atts::Dict{Any,Any},t::Union{DataType,Integer}=Float64)=NcVar{typeof(t)==DataType ? t : nctype2jltype[t],length(dimin)}(-1,-1,length(dimin),length(atts), typeof(t)==DataType ? jltype2nctype[t] : t,name,Array(Int,length(dimin)),dimin,atts,-1)

#Array methods
@generated function Base.size{T,N}(a::NcVar{T,N})
  :(@ntuple($N,i->Int(a.dim[i].dimlen)))
end
typealias IndR Union{Integer,UnitRange,Colon}
typealias ArNum Union{AbstractArray,Number}

Base.linearindexing(::NcVar)=Base.LinearSlow()
Base.getindex{T}(v::NcVar{T,1},i1::IndR)=readvar(v,i1)
Base.getindex{T}(v::NcVar{T,2},i1::IndR,i2::IndR)=readvar(v,i1,i2)
Base.getindex{T}(v::NcVar{T,3},i1::IndR,i2::IndR,i3::IndR)=readvar(v,i1,i2,i3)
Base.getindex{T}(v::NcVar{T,4},i1::IndR,i2::IndR,i3::IndR,i4::IndR)=readvar(v,i1,i2,i3,i4)
Base.getindex{T}(v::NcVar{T,5},i1::IndR,i2::IndR,i3::IndR,i4::IndR,i5::IndR)=readvar(v,i1,i2,i3,i4,i5)
Base.getindex{T}(v::NcVar{T,6},i1::IndR,i2::IndR,i3::IndR,i4::IndR,i5::IndR,i6::IndR)=readvar(v,i1,i2,i3,i4,i5,i6)

Base.setindex!{T}(v::NcVar{T,1},x::ArNum,i1::IndR)=putvar(v,x,i1)
Base.setindex!{T}(v::NcVar{T,2},x::ArNum,i1::IndR,i2::IndR)=putvar(v,x,i1,i2)
Base.setindex!{T}(v::NcVar{T,3},x::ArNum,i1::IndR,i2::IndR,i3::IndR)=putvar(v,x,i1,i2,i3)
Base.setindex!{T}(v::NcVar{T,4},x::ArNum,i1::IndR,i2::IndR,i3::IndR,i4::IndR)=putvar(v,x,i1,i2,i3,i4)
Base.setindex!{T}(v::NcVar{T,5},x::ArNum,i1::IndR,i2::IndR,i3::IndR,i4::IndR,i5::IndR)=putvar(v,x,i1,i2,i3,i4,i5)
Base.setindex!{T}(v::NcVar{T,6},x::ArNum,i1::IndR,i2::IndR,i3::IndR,i4::IndR,i5::IndR,i6::IndR)=putvar(v,x,i1,i2,i3,i4,i5,i6)

"""
The type `NcFile` represents a link to a NetCDF file. I can be retrieved by `NetCDF.open` or by `NetCDF.create`.
"""
type NcFile
  ncid::Int32
  nvar::Int32
  ndim::Int32
  ngatts::Int32
  vars::Dict{ASCIIString,NcVar}
  dim::Dict{ASCIIString,NcDim}
  gatts::Dict
  nunlimdimid::Int32
  name::ASCIIString
  omode::UInt16
  in_def_mode::Bool
end
#Define getindex method to retrieve a variable
Base.getindex(nc::NcFile,i::AbstractString)=haskey(nc.vars,i) ? nc.vars[i] : error("NetCDF file $(nc.name) does not have a variable named $(i)")

include("netcdf_helpers.jl")

const currentNcFiles=Dict{ASCIIString,NcFile}()

readvar!(nc::NcFile, varname::AbstractString, retvalsa::Array;start::Vector=defaultstart(nc[varname]),count::Vector=defaultcount(nc[varname]))=readvar!(nc[varname],retvalsa,start=start,count=count)

"""
readvar!(v::NcVar, d::Array;start::Vector=ones(UInt,ndims(d)),count::Vector=size(d))
is the mutating form of `readvar` which expects a pre-allocated array d, where the data are written to. It reads the values of the NcVar object `v`
writing the values to the preallocated array `d`.
If only parts of the variable are to be read, you can provide optionally `start` and `count`, which enable you to read blocks of data. `start` and `count` have the same
length as the number of variable dimensions. start gives the initial index for each dimension, while count gives the number of indices to be read along each dimension.
As a special case, setting a value in count to -1 will cause the function to read all values along this dimension.
"""
function readvar!(v::NcVar, retvalsa::Array;start::Vector=defaultstart(v),count::Vector=defaultcount(v))

  length(start) == v.ndim || error("Length of start ($(length(start))) must equal the number of variable dimensions ($(v.ndim))")
  length(count) == v.ndim || error("Length of start ($(length(count))) must equal the number of variable dimensions ($(v.ndim))")

  p=preparestartcount(start,count,v)

  length(retvalsa) != p && error("Size of output array ($(length(retvalsa))) does not equal number of elements to be read ($p)!")

  nc_get_vara_x!(v.ncid,v.varid,gstart,gcount,retvalsa)

  retvalsa
end
#readvar{T<:Integer}(nc::NcFile,varname::AbstractString,start::Array{T,1},count::Array{T,1})=readvar(nc,varname,start=start,count=count)

readvar(nc::NcFile, varname::AbstractString;start::Vector=defaultstart(nc[varname]),count::Vector=defaultcount(nc[varname]))=readvar(nc[varname],start=start,count=count)

"""
readvar{T,N}(v::NcVar{T,N};start::Vector=ones(Int,v.ndim),count::Vector=zeros(Int,v.ndim))
`readvar` reads the values of the NcVar object `v` and returns them.
If only parts of the variable are to be read, you can provide optionally `start` and `count`, which enable you to read blocks of data. `start` and `count` have the same
length as the number of variable dimensions. start gives the initial index for each dimension, while count gives the number of indices to be read along each dimension.
As a special case, setting a value in count to -1 will cause the function to read all values along this dimension.
"""
function readvar{T,N}(v::NcVar{T,N};start::Vector=defaultstart(v),count::Vector=defaultcount(v))
    retvalsa = Array(T,count...)
    readvar!(v, retvalsa, start=start, count=count)
    return retvalsa
end

"""
    readvar{T,N}(v::NcVar{T,N},I::IndR...)
Reading data from NetCDF file with array-style indexing. Integers and UnitRanges are valid indices for each dimension.
"""
function readvar{T,N}(v::NcVar{T,N},I::IndR...)
    count=ntuple(i->counti(I[i],v.dim[i].dimlen),length(I))
    retvalsa = zeros(T,count)
    readvar!(v, retvalsa, I...)
    return retvalsa
end

# Here are some functions for array-style indexing readvar
#For single indices
@generated function readvar{T,N}(v::NcVar{T,N},I::Integer...)
  N==length(I) || error("Dimension mismatch")
  quote
    checkbounds(v,I...)
    @nexprs $N i->gstart[v.ndim+1-i]=I[i]-1
    nc_get_var1_x(v.ncid,v.varid,gstart,T)::T
  end
end

firsti(i::Integer,l::Integer)=i-1
counti(i::Integer,l::Integer)=1
firsti(r::UnitRange,l::Integer)=start(r)-1
counti(r::UnitRange,l::Integer)=length(r)
firsti(r::Colon,l::Integer)=0
counti(r::Colon,l::Integer)=Int(l)
# For ranges
@generated function readvar!{T,N}(v::NcVar{T,N}, retvalsa::Array,I::IndR...)

  N==length(I) || error("Dimension mismatch")

  quote
    checkbounds(v,I...)
    @nexprs $N i->gstart[v.ndim+1-i]=firsti(I[i],v.dim[i].dimlen)
    @nexprs $N i->gcount[v.ndim+1-i]=counti(I[i],v.dim[i].dimlen)
    p=1
    @nexprs $N i->p=p*gcount[v.ndim+1-i]
    length(retvalsa) != p && error(string("Size of output array ($(length(retvalsa))) does not equal number of elements to be read (",p,")!"))
    nc_get_vara_x!(v.ncid,v.varid,gstart,gcount,retvalsa)
    retvalsa
  end
end

for (t,ending,arname) in funext
    fname = symbol("nc_get_vara_$ending")
    fname1 = symbol("nc_get_var1_$ending")
    arsym=symbol(arname)
    @eval nc_get_vara_x!(ncid::Integer,varid::Integer,start::Vector{UInt},count::Vector{UInt},retvalsa::Array{$t})=$fname(ncid,varid,start,count,retvalsa)
    @eval nc_get_var1_x(ncid::Integer,varid::Integer,start::Vector{UInt},::Type{$t})=begin $fname1(ncid,varid,start,$(arsym)); $(arsym)[1] end
end

function nc_get_vara_x!(ncid::Integer,varid::Integer,start::Vector{UInt},count::Vector{UInt},retvalsa::Array{AbstractString})
  retvalsa_c=Array(Ptr{UInt8},length(retvalsa))
  nc_get_vara_string(ncid,varid,start,count,retvalsa_c)
  for i=1:length(retvalsa)
    retvalsa[i]=bytestring(retvalsa_c[i])
  end
  nc_free_string(length(retvalsa_c),retvalsa_c)
  retvalsa
end

function nc_get_var1_x(ncid::Integer,varid::Integer,start::Vector{UInt},::Type{AbstractString})
  retvalsa_c=Array(Ptr{UInt8},1)
  nc_get_var1_string(ncid,varid,start,retvalsa_c)
  retval=bytestring(retvalsa_c[1])
  nc_free_string(1,retvalsa_c)
  retval
end

function putatt(ncid::Integer,varid::Integer,atts::Dict)
  for a in atts
    name=a[1]
    val=a[2]
    nc_put_att(ncid,varid,name,val)
  end
end

function putatt(nc::NcFile,varname::AbstractString,atts::Dict)
  varid = haskey(nc.vars,varname) ? nc.vars[varname].varid : NC_GLOBAL
  chdef=false
  if (!nc.in_def_mode)
    chdef=true
    nc_redef(nc.ncid)
  end
  putatt(nc.ncid,varid,atts)
  chdef ? nc_enddef(nc.ncid) : nothing
end

"""
    ncputatt(nc::String,varname::String,atts::Dict)
Writes name-value pairs of attributes given in `atts` to the file `nc`. If `varname` is a valid variable indentifier, the attribute(s) will be written to this variable.
Otherwise, global attributes will be generated. This function automatically overwrites existing attributes.
"""
function ncputatt(nc::AbstractString,varname::AbstractString,atts::Dict)
  nc = haskey(currentNcFiles,abspath(nc)) ? currentNcFiles[abspath(nc)] : open(nc,mode=NC_WRITE)
  if (nc.omode==NC_NOWRITE)
    fil=nc.name
    close(nc)
    println("reopening file in WRITE mode")
    open(fil,mode=NC_WRITE)
  end
  putatt(nc,varname,atts)
end

putvar(nc::NcFile,varname::AbstractString,vals::Array;start=ones(Int,length(size(vals))),count=[size(vals)...])=putvar(nc[varname],vals,start=start,count=count)

"""
    putvar(v::NcVar,vals::Array;start::Vector=ones(Int,length(size(vals))),count::Vector=[size(vals)...])
This function writes the values from the array `vals` to a netcdf file. `v` is the NcVar handle of the respective variable and `vals` an array
with the same dimension as the variable in the netcdf file. The optional parameter start gives the first index in each dimension along which the writing should
begin. It is assumed that the input array vals has the same number of dimensions as the and writing happens along these dimensions. However, you can specify
the number of values to be written along each dimension by adding an optional count argument, which is a vector whose length equals the number of variable dimensions.
"""
function putvar(v::NcVar,vals::Array;start::Vector=ones(Int,length(size(vals))),count::Vector=[size(vals)...])
  p=preparestartcount(start,count,v)
  nc_put_vara_x(v.ncid,v.varid,gstart,gcount,vals)
end

@generated function putvar{T,N}(v::NcVar{T,N},val::Any,I::IndR...)

  N==length(I) || error("Dimension mismatch")

  quote
    checkbounds(v,I...)
    @nexprs $N i->gstart[v.ndim+1-i]=firsti(I[i],v.dim[i].dimlen)
    @nexprs $N i->gcount[v.ndim+1-i]=counti(I[i],v.dim[i].dimlen)
    p=1
    @nexprs $N i->p=p*gcount[v.ndim+1-i]
    length(val) != p && error(string("Size of output array ($(length(retvalsa))) does not equal number of elements to be read (",p,")!"))
    nc_put_vara_x(v.ncid,v.varid,gstart,gcount,val)
  end
end

function putvar{T,N}(v::NcVar{T,N},val::Any,I::Integer...)

    N==length(I) || error("Dimension mismatch")
    quote
      checkbounds(v,I...)
      @nexprs $N i->gstart[v.ndim+1-i]=I[i]-1
      nc_put_var1_x(v.ncid,v.varid,gstart,val)
    end

end

for (t,ending,arname) in funext
    fname = symbol("nc_put_vara_$ending")
    fname1= symbol("nc_put_var1_$ending")
    arsym=symbol(arname)
    @eval nc_put_vara_x(ncid::Integer, varid::Integer, start, count, vals::Array{$t})=$fname(ncid,varid,start,count,vals)
    @eval nc_put_var1_x(ncid::Integer,varid::Integer,start::Vector{UInt},val::$t)=begin $(arsym)[1]=val; $fname1(ncid,varid,start,$(arsym)) end
end

function nc_put_vara_x{T<:AbstractString}(ncid::Integer, varid::Integer, start, count,vals::Array{T})
  vals_p=map(x->pointer(x.data),vals)
  nc_put_vara_string(ncid,varid,start,count,vals_p)
end

function nc_put_var1_x(ncid::Integer,varid::Integer,start::Vector{UInt},val::AbstractString)
  val_p=fill(pointer(val.data),1)
  nc_put_var1_string(ncid,varid,start,val_p)
end

"Synchronizes the changes made to the file and writes changes to the disk. If the argument is omitted, all open files are synchronized. "
function ncsync()
  for ncf in currentNcFiles
    _nc_sync_c(ncf[2].ncid)
  end
end

function sync(nc::NcFile)
  _nc_sync_c(nc.ncid)
end

#Function to close netcdf files
"Closes the file and writes changes to the disk. If argument is omitted, all open files are closed.   "
function ncclose(fil::AbstractString)
  if (haskey(currentNcFiles,abspath(fil)))
    close(currentNcFiles[abspath(fil)])
  else
    println("File $fil not currently opened.")
  end
end

function ncclose()
  for f in keys(currentNcFiles)
    ncclose(f)
  end
end

function setcompression(v::NcVar,mode)
    if v.compress > -1
        if (NC_NETCDF4 & mode)== 0
            warn("Compression only possible for NetCDF4 file format. Compression will be ingored.")
            v.compress=-1
        else
            v.compress=max(v.compress,9)
            nc_def_var_deflate(v.ncid,v.varid,Int32(1),Int32(1),v.compress);
        end
    end
end
"""
    NetCDF.create(name::String,varlist::Array{NcVar};gatts::Dict{Any,Any}=Dict{Any,Any}(),mode::UInt16=NC_NETCDF4)
This creates a new NetCDF file. Here, `filename`
 is the name of the file to be created and `varlist` an array of `NcVar` holding the variables that should appear in the file. In the optional
argument `gatts` you can specify a Dict containing global attributes and `mode` is the file type you want to create (NC_NETCDF4, NC_CLASSIC_MODEL or NC_64BIT_OFFSET).
"""
function create(name::AbstractString,varlist::Array{NcVar};gatts::Dict{Any,Any}=Dict{Any,Any}(),mode::UInt16=NC_NETCDF4)

  #Create the file
  id = nc_create(name,mode)

  # Collect Dimensions and set NetCDF ID
  vars=Dict{ASCIIString,NcVar}();
  dims=Set{NcDim}();
  for v in varlist
      v.ncid=id
      for d in v.dim
          push!(dims,d);
      end
  end
  nunlim=0
  ndim=Int32(length(dims))

  #Create the NcFile Object
  nc=NcFile(id,Int32(length(vars)),ndim,zero(Int32),vars,Dict{ASCIIString,NcDim}(),Dict{Any,Any}(),zero(Int32),name,NC_WRITE,true)

  for d in dims

      create_dim(nc,d)

    #Create dimension variable, if necessary
    if (length(d.vals)>0) & (!haskey(vars,d.name))
      push!(varlist,NcVar{Float64,1}(id,varida[1],1,length(d.atts),NC_DOUBLE,d.name,[d.dimid],[d],d.atts,-1))
    end

  end

  # Create variables in the file
  for v in varlist
      create_var(nc,v,mode)
  end

  # Put global attributes
  if !isempty(gatts)
    putatt(id,NC_GLOBAL,gatts)
  end

  # Leave define mode
  nc_enddef(nc)

  currentNcFiles[abspath(nc.name)]=nc

  for d in nc.dim
    #Write dimension variable
    if (length(d[2].vals)>0)
      putvar(nc,d[2].name,d[2].vals)
    end
  end

  return(nc)
end

create(name::AbstractString,varlist::NcVar...;gatts::Dict{Any,Any}=Dict{Any,Any}(),mode::UInt16=NC_NETCDF4)=create(name,NcVar[varlist[i] for i=1:length(varlist)];gatts=gatts,mode=mode)

"""
    NetCDF.close(nc::NcFile)
closes a NetCDF file handle
"""
function close(nco::NcFile)
  #Close file
  nc_close(nco.ncid)
  delete!(currentNcFiles,abspath(nco.name))
  return nco.ncid
end

"""
    NetCDF.open(fil::AbstractString,v::AbstractString; mode::Integer=NC_NOWRITE, readdimvar::Bool=false)
opens a NetCDF variable `v` in the NetCDF file `fil` and returns an `NcVar` handle that implements the AbstractArray interface for reading and writing.
"""
function open(fil::AbstractString,v::AbstractString; mode::Integer=NC_NOWRITE, readdimvar::Bool=false)
    nc=open(fil,mode=mode,readdimvar=readdimvar)
    nc.vars[v]
end

"""
    open(fil::AbstractString; mode::Integer=NC_NOWRITE, readdimvar::Bool=false)
opens the NetCDF file `fil` and returns a `NcFile` handle. The optional argument mode determines the `mode` in which the files is opened (`NC_NOWRITE` or `NC_WRITE`).
If you set `readdimvar=true`, then the dimension variables will be read when opening the file and added to the NcFIle object.
"""
function open(fil::AbstractString; mode::Integer=NC_NOWRITE, readdimvar::Bool=false)
  # Open netcdf file
  ncid=nc_open(fil,mode)

  #Get initial information
  ndim,nvar,ngatt,nunlimdimid = nc_inq(ncid)

  #Create ncdf object
  ncf=NcFile(ncid,Int32(nvar-ndim),ndim,ngatt,Dict{ASCIIString,NcVar}(),Dict{ASCIIString,NcDim}(),Dict{Any,Any}(),nunlimdimid,abspath(fil),mode,false)

  #Read global attributes
  ncf.gatts=getatts_all(ncid,NC_GLOBAL,ngatt)

  #Read dimensions
  for dimid = 0:ndim-1
    (name,dimlen)=nc_inq_dim(ncid,dimid)
    ncf.dim[name]=NcDim(ncid,dimid,-1,name,dimlen,[],Dict{Any,Any}())
  end

  #Read variable information
  for varid = 0:(nvar-1)
    (name,nctype,dimids,natts,vndim,isdimvar)=nc_inq_var(ncf,varid)
    if (isdimvar)
      ncf.dim[name].varid=varid
    end
    atts=getatts_all(ncid,varid,natts)
    vdim=Array(NcDim,length(dimids))
    i=1;
    for did in dimids
      vdim[i]=ncf.dim[getdimnamebyid(ncf,did)]
      i=i+1
    end
    ncf.vars[name]=NcVar{nctype2jltype[nctype],Int(vndim)}(ncid,Int32(varid),vndim,natts,nctype,name,dimids[vndim:-1:1],vdim[vndim:-1:1],atts,0)
  end
  readdimvar==true ? _readdimvars(ncf) : nothing
  currentNcFiles[abspath(ncf.name)]=ncf
  return ncf
end

# Define some high-level functions
"""
ncread(filename, varname, start=[1,1,...], count=[-1,-1,...] )
reads the values of the variable varname from file filename and returns the values in an array.
If only parts of the variable are to be read, you can provide optionally `start` and `count`, which enable you to read blocks of data. `start` and `count` have the same
length as the number of variable dimensions. start gives the initial index for each dimension, while count gives the number of indices to be read along each dimension.
As a special case, setting a value in count to -1 will cause the function to read all values along this dimension.
"""
function ncread{T<:Integer}(fil::AbstractString,vname::AbstractString;start::Array{T}=Array(Int,0),count::Array{T}=Array(Int,0))
  nc = haskey(currentNcFiles,abspath(fil)) ? currentNcFiles[abspath(fil)] : open(fil)
  length(start)==0 && (start=defaultstart(nc[vname]))
  length(count)==0 && (count=defaultcount(nc[vname]))
  x  = readvar(nc[vname],start=start,count=count)
  return x
end
ncread{T<:Integer}(fil::AbstractString,vname::AbstractString,start::Array{T,1},count::Array{T,1})=ncread(fil,vname,start=start,count=count)

"""
ncread! ( filename, varname, d, start=[1,1,...], count=[-1,-1,...] )
is the mutating form of `ncread` which expects a pre-allocated array d, where the data are written to. It reads the values of the variable varname from file filename
writing the values to the preallocated array `d`.
If only parts of the variable are to be read, you can provide optionally `start` and `count`, which enable you to read blocks of data. `start` and `count` have the same
length as the number of variable dimensions. start gives the initial index for each dimension, while count gives the number of indices to be read along each dimension.
As a special case, setting a value in count to -1 will cause the function to read all values along this dimension.
"""
function ncread!{T<:Integer}(fil::AbstractString,vname::AbstractString,vals::Array;start::Array{T}=ones(Int,ndims(vals)),count::Array{T}=Array(Int,size(vals)))
  nc = haskey(currentNcFiles,abspath(fil)) ? currentNcFiles[abspath(fil)] : open(fil)
  x  = readvar!(nc,vname,vals,start=start,count=count)
  return x
end

"prints information on the variables, dimension and attributes conatained in the file"
function ncinfo(fil::AbstractString)
  nc = haskey(currentNcFiles,abspath(fil)) ? currentNcFiles[abspath(fil)] : open(fil)
  return(nc)
end

iswritable(nc::NcFile)=(nc.omode & NC_WRITE) != zero(UInt16)

#High-level functions for writing data to a file
"""
    ncwrite{T<:Integer}(x::Array,fil::AbstractString,vname::AbstractString;start::Array{T,1}=ones(Int,length(size(x))),count::Array{T,1}=[size(x)...])
Writes the array `x` to the file `fil` and variable `vname`. If no `start` argument is supplied, writing starts at index 1 in each dimension.
You can supply the argument `start`, a vector that has the same number as the number of variable dimensions, that provides the indices where to start writing the data.
As default the number of values written along each dimension equals the dimension of the input array. However you can specify the along which dimension the data will be
written by setting a `count` argument, an integer vector indicating the number of values written along each dimension.
"""
function ncwrite{T<:Integer}(x::Array,fil::AbstractString,vname::AbstractString;start::Array{T,1}=ones(Int,length(size(x))),count::Array{T,1}=[size(x)...])
  nc = haskey(currentNcFiles,abspath(fil)) ? currentNcFiles[abspath(fil)] : open(fil,mode=NC_WRITE)
  if (nc.omode==NC_NOWRITE)
    close(nc)
    println("reopening file in WRITE mode")
    open(fil,mode=NC_WRITE)
  end
  putvar(nc,vname,x,start=start,count=count)
end
ncwrite(x::Array,fil::AbstractString,vname::AbstractString,start::Array)=ncwrite(x,fil,vname,start=start)

"This reads an attribute from the specified file and variable. To read global attributes, set varname to `Global`."
function ncgetatt(fil::AbstractString,vname::AbstractString,att::AbstractString)
  nc= haskey(currentNcFiles,abspath(fil)) ? currentNcFiles[abspath(fil)] : open(fil,NC_WRITE)
  return ( haskey(nc.vars,vname) ? get(nc.vars[vname].atts,att,nothing) : get(nc.gatts,att,nothing) )
end

#High-level function for creating files and variables
#
# if the file does not exist, it will be created
# if the file already exists, the variable will be added to the file

function create_dim(nc,dim)
    nc_redef(nc)
    nc_def_dim(nc.ncid,dim.name,dim.dimlen,dima);
    dim.dimid=dima[1];
    if !isempty(dim.vals)
      dumids[1]=dim.dimid
      nc_def_var(nc.ncid,dim.name,NC_DOUBLE,1,dumids,vara)
      nc.vars[dim.name]=NcVar{Float64,1}(nc.ncid,vara[1],1,0,NC_DOUBLE,dim.name,Int32[dim.dimid],NcDim[dim],Dict{Any,Any}(),-1)
    end
    nc.dim[dim.name]=dim;
end

function create_var(nc,v,mode)

    nc_redef(nc)

    v.dimids=Int32[v.dim[i].dimid for i=1:length(v.dim)]
    for i=1:v.ndim dumids[i]=v.dimids[v.ndim+1-i] end

    nc_def_var(nc.ncid,v.name,v.nctype,v.ndim,dumids,vara)

    v.varid=vara[1];
    nc.vars[v.name]=v;
    putatt(nc.ncid,v.varid,v.atts)
    setcompression(v,mode)
end

"""
    nccreate (filename, varname, dimensions ..., atts=atts,gatts=gatts,compress=compress,t=t,mode=mode)
This creates a variable in an existing netcdf file or creates a new file. `filename` and `varname` are strings.
After that follows a list of dimensions. Each dimension entry starts with a dimension name (a String), and
may be followed by a dimension length, an array with dimension values or a Dict containing dimension attributes.
Then the next dimension is entered and so on. Have a look at examples/high.jl for an example use.

Possible optional arguments are:

- **atts** Dict of attribute names and values to be assigned to the variable created
- **gatts** Dict of attribute names and values to be written as global attributes
- **compress** Integer [0..9] setting the compression level of the file, only valid if mode=NC_NETCDF4
- **t** variable type, currently supported types are: const NC_BYTE, NC_CHAR, NC_SHORT, NC_INT, NC_FLOAT, NC_LONG, NC_DOUBLE
- **mode** file creation mode, only valid when new file is created, choose one of: NC_NETCDF4, NC_CLASSIC_MODEL, NC_64BIT_OFFSET
"""
function nccreate(fil::AbstractString,varname::AbstractString,dims...;atts::Dict=Dict{Any,Any}(),gatts::Dict=Dict{Any,Any}(),compress::Integer=-1,t::Union{DataType,Integer}=NC_DOUBLE,mode::UInt16=NC_NETCDF4)
    # Checking dims argument for correctness
    dim=parsedimargs(dims)
    # open the file
    # create the NcVar object
    v=NcVar(varname,dim,atts=atts,compress=compress,t=t)
    # Test if the file already exists
    if (isfile(fil))
        nc=haskey(currentNcFiles,abspath(fil)) ? currentNcFiles[abspath(fil)] : open(fil,mode=NC_WRITE)
        if (nc.omode==NC_NOWRITE)
            close(nc)
            println("reopening file in WRITE mode")
            open(fil,NC_WRITE)
        end
        v.ncid = nc.ncid
        haskey(nc.vars,varname) ? error("Variable $varname already exists in file fil") : nothing
        # Check if dimensions exist, if not, create
        # Remember if dimension was created
        dcreate=Array(Bool,length(dim))
        for i=1:length(dim)
            if !haskey(nc.dim,dim[i].name)
                create_dim(nc,dim[i])
                v.dimids[i]=dim[i].dimid;
                dcreate[i] = true
            else
                v.dimids[i]=nc.dim[dim[i].name].dimid;
                v.dim[i] = nc.dim[dim[i].name]
                dcreate[i] = false
            end
        end
        # Create variable
        create_var(nc,v,mode)
        nc_enddef(nc)
        for i=1:length(dim)
            if dcreate[i] & !isempty(dim[i].vals)
                ncwrite(dim[i].vals,fil,dim[i].name)
            end
        end
    else
        nc = create(fil,v,gatts=gatts,mode=mode | NC_NOCLOBBER)
        for d in dim
            !isempty(d.vals) && ncwrite(d.vals,fil,d.name)
        end
    end
end

#show{T<:Any,N}(io::IO,a::NcVar{T,N})=println(io,a.name)
#showcompact{T<:Any,N}(io::IO,a::NcVar{T,N})=println(io,a.name)
function show(io::IO,nc::NcFile)
    nrow,ncol=Base.tty_size()
    hline = repeat("-",ncol)
    l1=div(ncol,3)
    l2=2*l1
    println(io,"")
    println(io,"##### NetCDF File #####")
    println(io,"")
    println(io,nc.name)
    println(io,"")
    println(io,"##### Dimensions #####")
    println(io,"")
    println(io,tolen("Name",l2),tolen("Length",l1))
    println(hline)
  for d in nc.dim
    println(io,tolen(d[2].name,l2),tolen(d[2].dimlen,l1))
  end
  l1=div(ncol,5)
  l2=2*l1
  println(io,"")
  println(io,"##### Variables #####")
  println(io,"")
  println(io,tolen("Name",l2),tolen("Type",l1),tolen("Dimensions",l2))
  println(hline)
  for v in nc.vars
    s1=string(tolen(v[2].name,l2))
    s2=string(tolen(nctype2string[Int(v[2].nctype)],l1))
    s3=""
    for d in v[2].dim
      s3=string(s3,d.name," ")
    end
    println(s1,s2,tolen(s3,l2))
  end
  l1=div(ncol,4)
  l2=2*l1
  println(io,"")
  println(io,"##### Attributes #####")
  println(io,"")
  println(io,tolen("Variable",l1),tolen("Name",l1),tolen("Value",l2))
  println(hline)
  for a in nc.gatts
    println(io,tolen("global",l1),tolen(a[1],l1),tolen(a[2],l2))
  end
  for v in nc.vars
    for a in v[2].atts
      println(io,tolen(v[2].name,l1),tolen(a[1],l1),tolen(a[2],l2))
    end
  end

end
tolen(s::Any,l::Number)=tolen(string(s),round(Int,l))
function tolen(s::AbstractString,l::Integer)
  ls = length(s)
  if ls<l
    return rpad(s,l)
  elseif ls>l
    return string(s[1:prevind(s,l-1)],"..")
  else
    return s
  end
end

end # Module


text:
module LPC

#
# TODO:
# [v] lpc-burg
# [ ] lattice-filter synth
# [ ] warped-lpc
# [ ] lpc-acorr
# [ ] lpc-covar
# [ ] Line-Spectral-Pairs
# [v] power-spectrum est
#

export lpc_burg, pow_spect
export log_area_ratio, inv_sine_coeffs

function lpc_burg{T <: Number}( x::AbstractVector{T}, p::Int)
# LPC (Linear-Predictive-Code) estimation, using the Burg-method
# @param: x - signal to predict
# @param: p - prediction order
# @retval: a - direct form coeffs
# @retval: r - reflection coefficients (usable for lattice-filter?)
# such that:
# minimizing x_prediction(n) = a*x
# ie: an IIR:
#  y = filt(1, a, white)
# should have similiar spectrum to x
# where 'white' is a signal with white spectrum
#
# The Burg-Method, which, provide both a minimal-phase (stable IIR)
# and estimates 'a', on a finite support [1]
# Combining the advantages on both traditional LPC methods
# (autocorrelation - minPhase, and covariance - finite-support)
#
# Implementing the loop from the article:
# [1] - ENHANCED PARTIAL TRACKING USING LINEAR PREDICTION
# (DAFX'03 article, Lagrange et al)
# http://www.elec.qmul.ac.uk/dafx03/proceedings/pdfs/dafx19.pdf
#
# Note: After implementing - I found matlab has an 'arburg()' - in
# the 'signal-processing-toolbox', with virtually identical code
# but hey, why use a library when you can easily reimplement it?
# anyhow, normalized interface.
#
# (c) Kobic, 2009, 2015, MIT License
# Translated my old(er) matlab code into julia::

# TODO: check memory allocations, gc etc!

    ef = x    # forward error
    eb = x    # backwards error
    a = [1; zeros(T, p)]  # prediction coeffs
    refl = zeros(T, p) # reflection coeffs
    k = zero(T)

# zero-mean - or not?
    prediction_err = dot(x,x) ./ length(x)  # variance

    for m in 1:p
        efp = ef[1:end-1]
        ebp = eb[2:end]
        k = -2 .* dot(ebp,efp) ./ ( dot(ebp,ebp) .+ dot(efp,efp) )
        refl[m] = k           # save reflection coeffs (if needed)
        ef = efp + k.*ebp
        eb = ebp + k.*efp
        a[1:m+1] = [a[1:m]; 0] + k.*[0; a[m:-1:1]]
        prediction_err *= (1 - k*k)
    end

    return a, prediction_err, refl
end

function pow_spect(ar, err, dft_order = 256)
    f = rfft( [ar; zeros(dft_order - length(ar))] )
    psd = err ./ (abs(f).^2 + eps() )

    return psd
end

function log_area_ratio(k)
    return log( (1 + k) ./ (1 - k) )
end

function inv_sine_coeffs(k)
    return 2asin(k) ./ pi
end

end # module


text:
# # Basics: Arithmetic, Assignment, and Strings
#
# Topics:
# 1. Infix arithmetic operators
# 2. Numeric literals
# 3. Comparisons
# 4. Assignments
# 5. Strings

#-

# ## Arithmetic should largely be familiar

3 + 7 # addition

#-

10 - 3 # subtraction

#-

20 * 5 # multiplication

#-

100 / 10 # division

#-

10 ^ 2 # exponentiation

#-

101 % 2 # remainder (modulus)

#-

sqrt(2) # square root

#-

√2 # Unicode to the rescue: \sqrt + TAB

# Note that dividing two integers yields a floating point number. There are two additional operators that may be helpful here:

10 / 6

#-

10 ÷ 6 # \div TAB or the `div` function

#-

div(10, 6)

#-

10 // 6

# ### Numbers: Many different ways to write the number forty-two

fortytwos = (42, 42.0, 4.20e1, 4.20f1, 84//2, 0x2a)

#-

for x in fortytwos
    show(x)
    println("\tisa $(typeof(x))")
end

# ### Bitwise arithmetic

0x2a & 0x70 # AND

#-

0x2a | 0x70 # OR

#-

42 & 112

#-

0b0010 << 2 # == 0b1000

# ### Logical operators

false && true # AND

#-

false || true # OR

# Note that they "short-circuit!"

x = -42
x > 0 || error("x must be positive")

# ### Comparisons

1 == 1.0 # Equality

#-

1 === 1.0 # Programmatically identical

#-

3 < π

#-

1 <= 1

#-

.1 + .2

#-

.1 + .2 ≈ .3 # \approx + TAB

# Comparisons "chain"
#
# Try inserting parentheses around one of these comparisons

2 == 2.0 == 0x02

#-

x = 42
0 < x < 100 || error("x must be between 0 and 100")

# ### Higher precision

big(2)^1000

big(pi)

big"0.1"

big(0.1)

# # Assignment
#
# Assignment in Julia is done with the single `=`. All it does is associates a name (on the left) to a value (on the right).

x = 1 # Use the name `x` for the value `1`

#-

y = x # Use the name `y` for the value that `x` refers to

#-

x = "hello!" # Decide you have a better use for the name `x`

#-

y # Is still the value 1

#-

# "Simultaneous" multiple assignment

x, y = y, x  # swap x and y
x, y = y, x  # swap back

#-

ϵ = eps(1.0) # You can make your own unicode names

#-

5ϵ # Juxtaposition is multiplication

#-

# We make use of juxtaposition for complex numbers

(1 + 2im)^2

# ## Updating operators
#
# All the infix arithmetic operators above can be used as "updating" operators in conjunction with an assignment:

y += 1

#-

## This is exactly the same as:
y = y + 1

# Note that it's just re-purposing the _same name_ for a new value. This means that the type might even change:

y /= 2

# # Strings

s1 = "I am a string."

#-

s2 = """I am also a string. """

#-

"Here, we get an "error" because it's ambiguous where this string ends "

#-

"""Look, Mom, no "errors"!!! """

#-

    println("""The other nice thing about triple-quoted
               string literals is that they ignore leading
               indentation, which is nice for long strings
               in real code. Try changing these quotes!""")

# Strings are not written with single `'`s — that's used for a single character:

first(s1)

#-

'⊂'

#-

'If you try writing a string in single-quotes, you will get an error'

# ## String interpolation

#-

# You can use the dollar sign inside a string to evaluate a Julia expression inside a string — either a single variable or a more complicated expression:

name = "Jane"
num_fingers = 10
num_toes = 10
println("Hello, my name is $name.")
println("I have $num_fingers fingers and $num_toes toes.")

#-

 println("That is $(num_fingers + num_toes) digits in all!!")

#-

@assert days == 365
@assert days_float == 365.0

max_stars_repo_path: Courses/Introduction to Julia/01. Arithmetic, Assignment, and Strings.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
module TestPoly
using Test

using Catlab.CategoricalAlgebra: is_isomorphic
using Poly

p = FinPoly([3,2,0])
@test p == FinPoly(FinFunction([1,1,1,2,2],3))
q = FinPoly([2,1])
@test p+q == FinPoly([3,2,0,2,1])
@test p * FinPoly([0]) == p
@test is_isomorphic(p * FinPoly([1]), FinPoly([4,3,1]))
@test is_isomorphic(p ⊗ FinPoly([1,2]), FinPoly([3,2,0,6,4,0]))

end


text:
# # 220: 2D Nonlinear Poisson with boundary reaction and boundary species
# ([source code](SOURCE_URL))

module Example220_NonlinearPoisson2D_BoundarySpecies

using Printf
using VoronoiFVM
using ExtendableGrids
using GridVisualize

function main(;n=10,Plotter=nothing,verbose=false,unknown_storage=:sparse)

    h=1.0/convert(Float64,n)
    X=collect(0.0:h:1.0)
    Y=collect(0.0:h:1.0)

    grid=VoronoiFVM.Grid(X,Y)

    k=1.0
    eps=1.0
    physics=VoronoiFVM.Physics(
    num_species=3,
    breaction=function(f,u,node)
        if  node.region==2
            f[1]=k*(u[1]-u[3])
            f[3]=k*(u[3]-u[1])+ k*(u[3]-u[2])
            f[2]=k*(u[2]-u[3])
        end
    end,

    bstorage=function(f,u,node)
        if  node.region==2
            f[3]=u[3]
        end
    end,

    flux=function(f,u0,edge)
        u=unknowns(edge,u0)
        f[1]=eps*(u[1,1]-u[1,2])
        f[2]=eps*(u[2,1]-u[2,2])
    end,

    source=function(f,node)
        x1=node[1]-0.5
        x2=node[2]-0.5
        f[1]=exp(-20.0*(x1^2+x2^2))
    end,

    storage=function(f,u,node)
        f[1]=u[1]
        f[2]=u[2]
    end
    )

    sys=VoronoiFVM.System(grid,physics,unknown_storage=unknown_storage)

    enable_species!(sys,1,[1])
    enable_species!(sys,2,[1])
    enable_boundary_species!(sys,3,[2])

    function tran32!(a,b)
        a[1]=b[2]
    end

    bgrid2=subgrid(grid,[2],boundary=true,transform=tran32!)

    inival=unknowns(sys)
    inival.=0.0
    U=unknowns(sys)

    eps=1.0e-2

    control=VoronoiFVM.NewtonControl()
    control.verbose=verbose
    control.tol_linear=1.0e-5
    control.tol_relative=1.0e-5
    control.max_lureuse=0
    tstep=0.01
    time=0.0
    istep=0
    u5=0
    p=GridVisualizer(Plotter=Plotter,layout=(3,1))
    while time<1
        time=time+tstep
        solve!(U,inival,sys,control=control,tstep=tstep)
        inival.=U
        if verbose
            @printf("time=%g\n",time)
        end
        tstep*=1.0
        istep=istep+1
        U_bound=view(U[3,:],bgrid2)
        u5=U_bound[5]
        scalarplot!(p[1,1],grid,U[1,:],clear=true)
        scalarplot!(p[2,1],grid,U[2,:])
        scalarplot!(p[3,1],bgrid2,U_bound,show=true,flimits=(0,0.0025))
    end
    return u5
end

function test()
    main(unknown_storage=:sparse) ≈ 0.0020781361856598
    main(unknown_storage=:dense) ≈ 0.0020781361856598
end
end


text:
export NumberTheory

module NumberTheory

export Z_star, Eulerφ, continued_fraction, mod_inverse, rand_primeto, factor_a_power_b
export is_order, order_from_float, find_order

"""
    Z_star(N::Int) -> Vector

returns the Z* group elements of `N`, i.e. {x | gcd(x, N) == 1}
"""
Z_star(N::Int) = filter(i->gcd(i, N)==1, 0:N-1)
Eulerφ(N) = length(Z_star(N))

"""
    continued_fraction(ϕ, niter::Int) -> Rational

obtain `s` and `r` from `ϕ` that satisfies `|s/r - ϕ| ≦ 1/2r²`
"""
continued_fraction(ϕ, niter::Int) = niter==0 ? floor(Int, ϕ) : floor(Int, ϕ) + 1//continued_fraction(1/mod(ϕ, 1), niter-1)
continued_fraction(ϕ::Rational, niter::Int) = niter==0 || ϕ.den==1 ? floor(Int, ϕ) : floor(Int, ϕ) + 1//continued_fraction(1/mod(ϕ, 1), niter-1)

"""
    mod_inverse(x::Int, N::Int) -> Int

Return `y` that `(x*y)%N == 1`, notice the `(x*y)%N` operations in Z* forms a group and this is the definition of inverse.
"""
function mod_inverse(x::Int, N::Int)
    for i=1:N
        (x*i)%N == 1 && return i
    end
    throw(ArgumentError("Can not find the inverse, $x is probably not in Z*($N)!"))
end

"""
    is_order(r, x, N) -> Bool

Returns true if `r` is the order of `x`, i.e. `r` satisfies `x^r % N == 1`.
"""
is_order(r, x, N) = powermod(x, r, N) == 1

"""
    find_order(x::Int, N::Int) -> Int

Find the order of `x` by brute force search.
"""
function find_order(x::Int, N::Int)
    findfirst(r->is_order(r, x, N), 1:N)
end

"""
    rand_primeto(N::Int) -> Int

Returns a random number `2 ≦ x < N` that is prime to `N`.
"""
function rand_primeto(N::Int)
    while true
        x = rand(2:N-1)
        d = gcd(x, N)
        if d == 1
            return x
        end
    end
end

"""
    order_from_float(ϕ, x, L) -> Int

Estimate the order of `x` to `L`, `r`, from a floating point number `ϕ ∼ s/r` using the continued fraction method.
"""
function order_from_float(ϕ, x, L)
    k = 1
    rnum = continued_fraction(ϕ, k)
    while rnum.den < L
        r = rnum.den
        if is_order(r, x, L)
            return r
        end
        k += 1
        rnum = continued_fraction(ϕ, k)
    end
    return nothing
end

"""
    factor_a_power_b(N::Int) -> (Int, Int) or nothing

Factorize `N` into the power form `a^b`.
"""
function factor_a_power_b(N::Int)
    y = log2(N)
    for b = 2:ceil(Int, y)
        x = 2^(y/b)
        u1 = floor(Int, x)
        u1^b == N && return (u1, b)
        (u1+1)^b == N && return (u1+1, b)
    end
end

end


text:
using Niancat.Formatters
import Niancat.Formatters: format

format(::SlackFormatter, r::Correct) = "Ordet $(r.word.w) är korrekt!"
format(::SlackFormatter, r::PuzzleIs) = r.puzzle
format(::SlackFormatter, r::NoPuzzleSet) = "Det finns inget pussel än."
format(::SlackFormatter, r::Incorrect) = "Ordet $(r.word.w) finns inte med i SAOL."
format(::SlackFormatter, r::Rejected) = "Inte ett giltigt pussel."
format(::SlackFormatter, r::NewPuzzle) = "Nya pusslet är $(r.puzzle)."

format(::SlackFormatter, r::CorrectNotification) = "$(r.user) löste nian!"


text:
function attemptLoadConfig(fn::String, bufferTime::Number=0, default::Dict{K, V}=Dict{Any, Any}()) where {K, V}
    try
        return loadConfig(fn, bufferTime, default)

    catch
        # Doesn't really matter what it is failing here, we would still want the default config

        failedFn = fn * ".failed"
        warningMessage = """Failed to parse config at "$fn".
            Program will exit if defaults are not used.
            Use defaults and continue?"""

        if ask_dialog(warningMessage)
            mv(fn, failedFn, force=true)

            config = Config(fn, bufferTime, default)
            saveConfig(config)

            return config

        else
            exit()
        end
    end
end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
product_manifold_split(::Type{PM{M1,M2}}) where {M1,M2} = M1, M2

dim(::Union{PM{M1,M2},Type{PM{M1,M2}}}) where {M1,M2} = dim(M1) + dim(M2)
embdim(::Union{PM{M1,M2},Type{PM{M1,M2}}}) where {M1,M2} = embdim(M1) + embdim(M2)

dims(::Union{PM{M1,M2},Type{PM{M1,M2}}}) where {M1,M2} = dim(M1), dim(M2)
embdims(::Union{PM{M1,M2},Type{PM{M1,M2}}}) where {M1,M2} = embdim(M1), embdim(M2)

product_chart_split(C::Chart{Tuple{I,J},PM{M1,M2}}) where {M1,M2,I,J} =
    Chart{I,M1}(), Chart{J,M2}()
product_chart_merge(C1::Chart{I,M1}, C2::Chart{J,M2}) where {M1,M2,I,J} =
    Chart{Tuple{I,J},PM{M1,M2}}()

isglobal(::Chart{Tuple{I,J},PM{M1,M2}}) where {M1,M2,I,J} =
    isglobal(Chart{I,M1}()) && isglobal(Chart{J,M2}())

function product_manifold_emb_splitview(pe, ::Type{M}) where M <: ProductManifold
    m = embdim(M)
    m1, m2 = embdims(M)
    x1e, x2e = pe[SVector{m1}(1:m1)], pe[SVector{m2}(m1+1:m)]
end

function product_manifold_chart_splitview(p, ::Type{M}) where M <: ProductManifold
    m = dim(M)
    m1, m2 = dims(M)
    x1, x2 = p[SVector{m1}(1:m1)], p[SVector{m2}(m1+1:m)]
end

function emb_to_chart(xe, C::Chart{I,M}) where {M<:PM,I}
    C1, C2 = product_chart_split(C)
    x1e, x2e = product_manifold_emb_splitview(xe, M)
    x1, x2 = emb_to_chart(x1e, C1), emb_to_chart(x2e, C2)
    x = [x1; x2]
end

function chart_to_emb(x, C::Chart{I,M}) where {M<:PM,I}
    C1, C2 = product_chart_split(C)
    x1, x2 = product_manifold_chart_splitview(x, M)
    x1e, x2e = chart_to_emb(x1, C1), chart_to_emb(x2, C2)
    xe = [x1e; x2e]
end

function chart_transition(xa, Ca::Chart{Tuple{Ia,Ja},M},
        Cb::Chart{Tuple{Ib,Jb},M}) where {M<:PM,Ia,Ja,Ib,Jb}
    C1a, C2a = product_chart_split(Ca)
    C1b, C2b = product_chart_split(Cb)
    x1a, x2a = product_manifold_chart_splitview(xa, M)
    x1b, x2b = chart_transition(x1a, C1a, C1b), chart_transition(x2a, C2a, C2b)
    xb = [x1b; x2b]
end

function choose_chart_emb(xe, C::Chart{Tuple{I,J},M}) where {M<:PM,I,J}
    C1, C2 = product_chart_split(C)
    x1e, x2e = product_manifold_emb_splitview(xe, M)
    C1c, C2c = choose_chart_emb(x1e, C1), choose_chart_emb(x2e, C2)
    product_chart_merge(C1c, C2c)
end

function choose_chart_chart(x, C::Chart{Tuple{I,J},M}) where {M<:PM,I,J}
    C1, C2 = product_chart_split(C)
    x1, x2 = product_manifold_chart_splitview(x, M)
    C1c, C2c = choose_chart_chart(x1, C1), choose_chart_chart(x2, C2)
    product_chart_merge(C1c, C2c)
end


text:
julia> # Load DataFrame package

julia> require("DataFrames")

julia> using DataFrames

julia>

julia> # Load a CSV file into a DataFrame

julia> df = read_table(Pkg.dir("DataFrames", "demo", "toy_example.csv"))
6x3 DataFrame:
          A       B   C
[1,]    2.5   "One" 3.0
[2,]    3.6   "One" 5.0
[3,]    3.5   "Two" 3.0
[4,]    4.5   "Two" 5.0
[5,]    4.5 "Three" 3.0
[6,]    5.5 "Three" 5.0

julia>

julia> # Look at the data structure

julia> # dump() is like R's str()

julia> dump(df)
DataFrame  6 observations of 3 variables
  A: DataArray{Float64,1}(6) [2.5,3.6,3.5,4.5]
  B: DataArray{UTF8String,1}(6) ["One","One","Two","Two"]
  C: DataArray{Float64,1}(6) [3.0,5.0,3.0,5.0]

julia>

julia> # Look at the internal data structure

julia> idump(df)
DataFrame
  columns: Array(Any,(3,))
    1: DataArray{Float64,1}
      data: Array(Float64,(6,)) [2.5, 3.6, 3.5, 4.5, 4.5, 5.5]
      na: BitArray{1}
        chunks: Array(Uint64,(1,)) [0x0000000000000000]
        dims: Array(Int64,(1,)) [6]
    2: DataArray{UTF8String,1}
      data: Array(UTF8String,(6,)) ["One", "One", "Two", "Two", "Three", "Three"]
      na: BitArray{1}
        chunks: Array(Uint64,(1,)) [0x0000000000000000]
        dims: Array(Int64,(1,)) [6]
    3: DataArray{Float64,1}
      data: Array(Float64,(6,)) [3.0, 5.0, 3.0, 5.0, 3.0, 5.0]
      na: BitArray{1}
        chunks: Array(Uint64,(1,)) [0x0000000000000000]
        dims: Array(Int64,(1,)) [6]
  colindex: Index
    lookup: Dict{Union(ASCIIString,UTF8String),Union(AbstractArray{Real,1},Real)}
      slots: Array(Uint8,(16,)) [0x00, 0x00, 0x00, 0x00, 0x00, 0x00  …  0x00, 0x01, 0x00, 0x00, 0x00, 0x00]
      keys: Array(Union(ASCIIString,UTF8String),(16,)) [#undef, #undef, #undef, #undef  …  "C", #undef, #undef, #undef, #undef]
      vals: Array(Union(AbstractArray{Real,1},Real),(16,)) [#undef, #undef, #undef, #undef  …  3, #undef, #undef, #undef, #undef]
      ndel: Int64 0
      count: Int64 3
      deleter: identity
    names: Array(Union(ASCIIString,UTF8String),(3,)) ["A", "B", "C"]

julia>

julia> # Print out a summary of each column

julia> describe(df)
A
Min      2.5
1st Qu.  3.525
Median   4.05
Mean     4.016666666666667
3rd Qu.  4.5
Max      5.5

B
Length: 6
Type  : UTF8String
NAs   : 0

C
Min      3.0
1st Qu.  3.0
Median   4.0
Mean     4.0
3rd Qu.  5.0
Max      5.0

julia>

julia> # head/tail of a DataFrame

julia> head(df)
6x3 DataFrame:
          A       B   C
[1,]    2.5   "One" 3.0
[2,]    3.6   "One" 5.0
[3,]    3.5   "Two" 3.0
[4,]    4.5   "Two" 5.0
[5,]    4.5 "Three" 3.0
[6,]    5.5 "Three" 5.0

julia> head(df, 3)
3x3 DataFrame:
          A     B   C
[1,]    2.5 "One" 3.0
[2,]    3.6 "One" 5.0
[3,]    3.5 "Two" 3.0

julia> tail(df, 2)
2x3 DataFrame:
          A       B   C
[1,]    4.5 "Three" 3.0
[2,]    5.5 "Three" 5.0

julia>

julia> # Select all rows where column A is greater than 4.0

julia> # Element-wise operators in Julia usually have a leading "."

julia> df[:( A .> 4.0 ), :]
3x3 DataFrame:
          A       B   C
[1,]    4.5   "Two" 5.0
[2,]    4.5 "Three" 3.0
[3,]    5.5 "Three" 5.0

julia>

julia> # Make a new column using within

julia> df2 = within(df, :( D = A + C ))
6x4 DataFrame:
          A       B   C    D
[1,]    2.5   "One" 3.0  5.5
[2,]    3.6   "One" 5.0  8.6
[3,]    3.5   "Two" 3.0  6.5
[4,]    4.5   "Two" 5.0  9.5
[5,]    4.5 "Three" 3.0  7.5
[6,]    5.5 "Three" 5.0 10.5

julia>

julia> # This is similar, but now changes apply directly to df

julia> within!(df, quote
           D = A + C
           E = A + sum(C)
       end)
6x5 DataFrame:
          A       B   C    D    E
[1,]    2.5   "One" 3.0  5.5 26.5
[2,]    3.6   "One" 5.0  8.6 27.6
[3,]    3.5   "Two" 3.0  6.5 27.5
[4,]    4.5   "Two" 5.0  9.5 28.5
[5,]    4.5 "Three" 3.0  7.5 28.5
[6,]    5.5 "Three" 5.0 10.5 29.5

julia>

julia> dump(df)
DataFrame  6 observations of 5 variables
  A: DataArray{Float64,1}(6) [2.5,3.6,3.5,4.5]
  B: DataArray{UTF8String,1}(6) ["One","One","Two","Two"]
  C: DataArray{Float64,1}(6) [3.0,5.0,3.0,5.0]
  D: DataArray{Float64,1}(6) [5.5,8.6,6.5,9.5]
  E: DataArray{Float64,1}(6) [26.5,27.6,27.5,28.5]

julia>

julia> # Create a new DataFrame based on operations on another DataFrame

julia> # This is similar to plyr's summarise()

julia> df3 = based_on(df, quote
           ct = cut(replaceNA(A, 0.0), 3)
           sum_A = sum(A)
       end)
6x2 DataFrame:
                     ct sum_A
[1,]    "[2.5,3.56667]"  24.1
[2,]    "(3.56667,4.5]"  24.1
[3,]    "[2.5,3.56667]"  24.1
[4,]    "(3.56667,4.5]"  24.1
[5,]    "(3.56667,4.5]"  24.1
[6,]        "(4.5,5.5]"  24.1

julia>

julia> # cut makes a PooledDataVector that is like R's factors, but

julia> # PooledDataVectors can contain more than just strings. Here's

julia> # the internal structure of a PooledDataVector:

julia> idump(df3["ct"])
PooledDataVector{ASCIIString}
  refs: Array(Uint16,(6,)) [0x0001, 0x0002, 0x0001, 0x0002, 0x0002, 0x0003]
  pool: Array(ASCIIString,(3,)) ["[2.5,3.56667]", "(3.56667,4.5]", "(4.5,5.5]"]

julia>

julia> # In DataFrame, copies of data are minimized, especially for column

julia> # operations

julia> # These are both the same entity; change one, and you change the other:

julia> df2 = df
6x5 DataFrame:
          A       B   C    D    E
[1,]    2.5   "One" 3.0  5.5 26.5
[2,]    3.6   "One" 5.0  8.6 27.6
[3,]    3.5   "Two" 3.0  6.5 27.5
[4,]    4.5   "Two" 5.0  9.5 28.5
[5,]    4.5 "Three" 3.0  7.5 28.5
[6,]    5.5 "Three" 5.0 10.5 29.5

julia>

julia> colA = df2["A"]
6-element Float64 DataArray
 2.5
 3.6
 3.5
 4.5
 4.5
 5.5

julia> colA[1] = 99.0   # This changes df and df2, too.
99.0

julia> df
6x5 DataFrame:
           A       B   C    D    E
[1,]    99.0   "One" 3.0  5.5 26.5
[2,]     3.6   "One" 5.0  8.6 27.6
[3,]     3.5   "Two" 3.0  6.5 27.5
[4,]     4.5   "Two" 5.0  9.5 28.5
[5,]     4.5 "Three" 3.0  7.5 28.5
[6,]     5.5 "Three" 5.0 10.5 29.5

julia>


text:
function linspace(start::T, stop::T, np::I) where {T <: AbstractFloat,I <: Integer}
    d = (stop - start) / (np - 1)
    p = Vector{T}(undef, np)
    p[1] = start
    p[end] = stop
    for i = 2:np - 1
        p[i] = start + (i - 1) * d
    end
    return p
end

function cosspace(start::T, stop::T, np::I) where {T <: AbstractFloat,I <: Integer}
    p = linspace(zero(T), convert(T, π), np)
    r = (stop - start) / 2
    for i = 1:np
        p[i] = start + r * (one(T) - cos(p[i]))
    end
    return p
end

function gen_airfoil(a::AbstractAirfoil{T}, np::Int=65; cos_space::Bool=true, orient::Symbol=:CW) where {T <: AbstractFloat}
    if cos_space
        xc = cosspace(zero(T), one(T), np)
    else
        xc = linspace(zero(T), one(T), np)
    end
    xu = Vector{T}(undef, np)
    yu = Vector{T}(undef, np)
    xl = Vector{T}(undef, np)
    yl = Vector{T}(undef, np)
    @inbounds for i = 1:np
        xu[i] = x_upper(a, xc[i])
        yu[i] = y_upper(a, xc[i])
        xl[i] = x_lower(a, xc[i])
        yl[i] = y_lower(a, xc[i])
    end
    return if orient == :CW
        vcat(xl[end:-1:2],xu), vcat(yl[end:-1:2],yu)
    elseif orient == :CCW
        vcat(xu[end:-1:2],xl), vcat(yu[end:-1:2],yl)
    else
        throw(ArgumentError("invalid orient=$orient, must be :CW or :CCW"))
    end
end

function n_upper(a::AbstractAirfoil{T}, xc::T)::Tuple{T,T} where {T <: AbstractFloat}
    θ = atan(dy_upper(a, xc))
    return -sin(θ), cos(θ)
end

function n_lower(a::AbstractAirfoil{T}, xc::T)::Tuple{T,T} where {T <: AbstractFloat}
    θ = atan(dy_lower(a, xc))
    return sin(θ), -cos(θ)
end

@inline n_upper(a::AbstractAirfoil{T}, xc::R) where {T <: AbstractFloat,R <: Real} = n_upper(a, convert(T, xc))
@inline n_lower(a::AbstractAirfoil{T}, xc::R) where {T <: AbstractFloat,R <: Real} = n_lower(a, convert(T, xc))

function t_upper(a::AbstractAirfoil{T}, xc::T)::Tuple{T,T} where {T <: AbstractFloat}
    dy = dy_upper(a, xc)
    if isinf(dy)
        return zero(T), sign(dy)
    else
        n = sqrt(dy^2 + oneunit(T))
        return 1 / n, dy / n
    end
end

function t_lower(a::AbstractAirfoil{T}, xc::T)::Tuple{T,T} where {T <: AbstractFloat}
    dy = dy_lower(a, xc)
    if isinf(dy)
        return zero(T), sign(dy)
    else
        n = sqrt(dy^2 + oneunit(T))
        return 1 / n, dy / n
    end
end

@inline t_upper(a::AbstractAirfoil{T}, xc::R) where {T <: AbstractFloat,R <: Real} = t_upper(a, convert(T, xc))
@inline t_lower(a::AbstractAirfoil{T}, xc::R) where {T <: AbstractFloat,R <: Real} = t_lower(a, convert(T, xc))

for fname in [:x_lower, :x_upper, :y_upper, :y_lower, :dy_upper, :dy_lower]
    @eval begin
        @inline function ($fname)(af::AbstractAirfoil, xc::AbstractVector)
            y = similar(xc)
            @inbounds for i in eachindex(xc)
                y[i] = ($fname)(af, xc[i])
            end
            y
        end
    end
end

for fname in [:fy_upper, :fy_lower, :n_upper, :n_lower, :t_upper, :t_lower]
    @eval begin
        @inline function ($fname)(af::AbstractAirfoil, xc::AbstractVector)
            y = similar(xc, 2, length(xc))
            @inbounds for i in eachindex(xc)
                y[:,i] .= ($fname)(af, xc[i])
            end
            y
        end
    end
end

function normalize_airfoil_data(xu, yu, xl, yl)
    iu = argmin(xu)
    il = argmin(xl)
    @assert xu[1] == xl[1]
    @assert yu[1] == yl[1]
    @assert iu == 1 || il == 1
    if iu > 1
        xl1 = vcat(xu[iu:-1:2], xl)
        yl1 = vcat(yu[iu:-1:2], yl)
        xu1 = xu[iu:end]
        yu1 = yu[iu:end]
    elseif il > 1
        xu1 = vcat(xl[il:-1:2], xu)
        yu1 = vcat(yl[il:-1:2], yu)
        xl1 = xl[il:end]
        yl1 = yl[il:end]
    else
        xu1 = copy(xu)
        yu1 = copy(yu)
        xl1 = copy(xl)
        yl1 = copy(yl)
    end
    x0 = xl1[1]
    y0 = yl1[1]
    @. xu1 -= x0
    @. yu1 -= y0
    @. xl1 -= x0
    @. yl1 -= y0
    x1 = (xl1[end] + xu1[end]) / 2
    @. xu1 /= x1
    @. yu1 /= x1
    @. xl1 /= x1
    @. yl1 /= x1
    xu1, yu1, xl1, yl1
end


text:
using QuantumOptics

b = SpinBasis(1//2)

sx = sigmax(b)
sy = sigmay(b)
sz = sigmaz(b)

N = 5

spinchain = b ⊗ b

for ispin in 3:N
    global spinchain
    spinchain = spinchain ⊗ b
end

XX = sx ⊗ sx
YY = sy ⊗ sy
ZZ = sz ⊗ sz

Jz=1

HIsing = embed(spinchain, [1,2] , Jz*ZZ)
for ispin in 2:N-1
    global HIsing
    HIsing += embed(spinchain, [ispin,ispin+1] , Jz*ZZ)
end

# make sure H is numerically Hermitian
HIsing = ( HIsing + dagger(HIsing) )/2

# find first 5 eigenvalues using sparse matrix
E5 = eigenenergies(HIsing, 5)
println(" ")
println("Lowest five eigenvalues")
println(E5)
println(" ")

# full spectrum
EE, UU = eigenstates(dense(HIsing))
println(" ")
println("All Eigenvalues")
println(EE)
println(" ")

max_stars_repo_name: BlackBodyRadiation/JupyterTNN

max_issues_repo_name: BlackBodyRadiation/JupyterTNN

max_forks_repo_name: BlackBodyRadiation/JupyterTNN


text:
# # Langmuir turbulence example
#
# This example implements a Langmuir turbulence simulation reported in section
# 4 of
#
# > [Wagner et al., "Near-inertial waves and turbulence driven by the growth of swell", Journal of Physical Oceanography (2021)](https://journals.ametsoc.org/view/journals/phoc/51/5/JPO-D-20-0178.1.xml)
#
# This example demonstrates
#
#   * How to run large eddy simulations with surface wave effects via the
#     Craik-Leibovich approximation.
#
#   * How to specify time- and horizontally-averaged output.

# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add Oceananigans, CairoMakie"
# ```

using Oceananigans
using Oceananigans.Units: minute, minutes, hours

# ## Model set-up
#
# To build the model, we specify the grid, Stokes drift, boundary conditions, and
# Coriolis parameter.
#
# ### Domain and numerical grid specification
#
# We use a modest resolution and the same total extent as Wagner et al. 2021,

grid = RectilinearGrid(size=(32, 32, 32), extent=(128, 128, 64))

# ### The Stokes Drift profile
#
# The surface wave Stokes drift profile prescribed in Wagner et al. 2021,
# corresponds to a 'monochromatic' (that is, single-frequency) wave field.
#
# A monochromatic wave field is characterized by its wavelength and amplitude
# (half the distance from wave crest to wave trough), which determine the wave
# frequency and the vertical scale of the Stokes drift profile.

using Oceananigans.BuoyancyModels: g_Earth

 amplitude = 0.8 # m
wavelength = 60 # m
wavenumber = 2π / wavelength # m⁻¹
 frequency = sqrt(g_Earth * wavenumber) # s⁻¹

## The vertical scale over which the Stokes drift of a monochromatic surface wave
## decays away from the surface is `1/2wavenumber`, or
const vertical_scale = wavelength / 4π

## Stokes drift velocity at the surface
const Uˢ = amplitude^2 * wavenumber * frequency # m s⁻¹

# The `const` declarations ensure that Stokes drift functions compile on the GPU.
# To run this example on the GPU, write `architecture = GPU()` in the constructor
# for `NonhydrostaticModel` below.
#
# The Stokes drift profile is

uˢ(z) = Uˢ * exp(z / vertical_scale)

# which we'll need for the initial condition.
#
# !!! info "The Craik-Leibovich equations in Oceananigans"
#     Oceananigans implements the Craik-Leibovich approximation for surface wave effects
#     using the _Lagrangian-mean_ velocity field as its prognostic momentum variable.
#     In other words, `model.velocities.u` is the Lagrangian-mean ``x``-velocity beneath surface
#     waves. This differs from models that use the _Eulerian-mean_ velocity field
#     as a prognostic variable, but has the advantage that ``u`` accounts for the total advection
#     of tracers and momentum, and that ``u = v = w = 0`` is a steady solution even when Coriolis
#     forces are present. See the
#     [physics documentation](https://clima.github.io/OceananigansDocumentation/stable/physics/surface_gravity_waves/)
#     for more information.
#
# The vertical derivative of the Stokes drift is

∂z_uˢ(z, t) = 1 / vertical_scale * Uˢ * exp(z / vertical_scale)

# Finally, we note that the time-derivative of the Stokes drift must be provided
# if the Stokes drift and surface wave field undergoes _forced_ changes in time.
# In this example, the Stokes drift is constant
# and thus the time-derivative of the Stokes drift is 0.

# ### Boundary conditions
#
# At the surface at ``z=0``, Wagner et al. 2021 impose

Qᵘ = -3.72e-5 # m² s⁻², surface kinematic momentum flux

u_boundary_conditions = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ))

# Wagner et al. 2021 impose a linear buoyancy gradient `N²` at the bottom
# along with a weak, destabilizing flux of buoyancy at the surface to faciliate
# spin-up from rest.

Qᵇ = 2.307e-9 # m³ s⁻², surface buoyancy flux
N² = 1.936e-5 # s⁻², initial and bottom buoyancy gradient

b_boundary_conditions = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ),
                                                bottom = GradientBoundaryCondition(N²))

# !!! info "The flux convention in Oceananigans"
#     Note that Oceananigans uses "positive upward" conventions for all fluxes. In consequence,
#     a negative flux at the surface drives positive velocities, and a positive flux of
#     buoyancy drives cooling.

# ### Coriolis parameter
#
# Wagner et al. (2021) use

coriolis = FPlane(f=1e-4) # s⁻¹

# which is typical for mid-latitudes on Earth.

# ## Model instantiation
#
# We are ready to build the model. We use a fifth-order Weighted Essentially
# Non-Oscillatory (WENO) advection scheme and the `AnisotropicMinimumDissipation`
# model for large eddy simulation. Because our Stokes drift does not vary in ``x, y``,
# we use `UniformStokesDrift`, which expects Stokes drift functions of ``z, t`` only.

model = NonhydrostaticModel(; grid,
                            advection = WENO5(),
                            timestepper = :RungeKutta3,
                            tracers = :b,
                            buoyancy = BuoyancyTracer(),
                            coriolis = coriolis,
                            closure = AnisotropicMinimumDissipation(),
                            stokes_drift = UniformStokesDrift(∂z_uˢ=∂z_uˢ),
                            boundary_conditions = (u=u_boundary_conditions, b=b_boundary_conditions))

# ## Initial conditions
#
# We make use of random noise concentrated in the upper 4 meters
# for buoyancy and velocity initial conditions,

Ξ(z) = randn() * exp(z / 4)
nothing # hide

# Our initial condition for buoyancy consists of a surface mixed layer 33 m deep,
# a deep linear stratification, plus noise,

initial_mixed_layer_depth = 33 # m
stratification(z) = z < - initial_mixed_layer_depth ? N² * z : N² * (-initial_mixed_layer_depth)

bᵢ(x, y, z) = stratification(z) + 1e-1 * Ξ(z) * N² * model.grid.Lz

# The simulation we reproduce from Wagner et al. (2021) is zero Lagrangian-mean velocity.
# This initial condition is consistent with a wavy, quiescent ocean suddenly impacted
# by winds. To this quiescent state we add noise scaled by the friction velocity to ``u`` and ``w``.

u★ = sqrt(abs(Qᵘ))
uᵢ(x, y, z) = u★ * 1e-1 * Ξ(z)
wᵢ(x, y, z) = u★ * 1e-1 * Ξ(z)

set!(model, u=uᵢ, w=wᵢ, b=bᵢ)

# ## Setting up the simulation

simulation = Simulation(model, Δt=45.0, stop_time=4hours)

# We use the `TimeStepWizard` for adaptive time-stepping
# with a Courant-Freidrichs-Lewy (CFL) number of 1.0,

wizard = TimeStepWizard(cfl=1.0, max_change=1.1, max_Δt=1minute)

simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))

# ### Nice progress messaging
#
# We define a function that prints a helpful message with
# maximum absolute value of ``u, v, w`` and the current wall clock time.

using Printf

function print_progress(simulation)
    u, v, w = simulation.model.velocities

    ## Print a progress message
    msg = @sprintf("i: %04d, t: %s, Δt: %s, umax = (%.1e, %.1e, %.1e) ms⁻¹, wall time: %s\n",
                   iteration(simulation),
                   prettytime(time(simulation)),
                   prettytime(simulation.Δt),
                   maximum(abs, u), maximum(abs, v), maximum(abs, w),
                   prettytime(simulation.run_wall_time))

    @info msg

    return nothing
end

simulation.callbacks[:progress] = Callback(print_progress, IterationInterval(10))

# ## Output
#
# ### A field writer
#
# We set up an output writer for the simulation that saves all velocity fields,
# tracer fields, and the subgrid turbulent diffusivity.

output_interval = 5minutes

fields_to_output = merge(model.velocities, model.tracers, (; νₑ=model.diffusivity_fields.νₑ))

simulation.output_writers[:fields] =
    JLD2OutputWriter(model, fields_to_output,
                     schedule = TimeInterval(output_interval),
                     filename = "langmuir_turbulence_fields.jld2",
                     overwrite_existing = true)

# ### An "averages" writer
#
# We also set up output of time- and horizontally-averaged velocity field and
# momentum fluxes,

u, v, w = model.velocities
b = model.tracers.b

 U = Average(u, dims=(1, 2))
 V = Average(v, dims=(1, 2))
 B = Average(b, dims=(1, 2))
wu = Average(w * u, dims=(1, 2))
wv = Average(w * v, dims=(1, 2))

simulation.output_writers[:averages] =
    JLD2OutputWriter(model, (; U, V, B, wu, wv),
                     schedule = AveragedTimeInterval(output_interval, window=2minutes),
                     filename = "langmuir_turbulence_averages.jld2",
                     overwrite_existing = true)

# ## Running the simulation
#
# This part is easy,

run!(simulation)

# # Making a neat movie
#
# We look at the results by loading data from file with FieldTimeSeries,
# and plotting vertical slices of ``u`` and ``w``, and a horizontal
# slice of ``w`` to look for Langmuir cells.

using CairoMakie

time_series = (;
     w = FieldTimeSeries("langmuir_turbulence_fields.jld2", "w"),
     u = FieldTimeSeries("langmuir_turbulence_fields.jld2", "u"),
     B = FieldTimeSeries("langmuir_turbulence_averages.jld2", "B"),
     U = FieldTimeSeries("langmuir_turbulence_averages.jld2", "U"),
     V = FieldTimeSeries("langmuir_turbulence_averages.jld2", "V"),
    wu = FieldTimeSeries("langmuir_turbulence_averages.jld2", "wu"),
    wv = FieldTimeSeries("langmuir_turbulence_averages.jld2", "wv"))

times = time_series.w.times
xw, yw, zw = nodes(time_series.w)
xu, yu, zu = nodes(time_series.u)
nothing # hide

# Finally, we're ready to animate using Makie.
n = Observable(1)

wxy_title = @lift @sprintf("w(x, y, t) at z=-8 m and t = %s ", prettytime(times[$n]))
wxz_title = @lift @sprintf("w(x, z, t) at y=0 m and t = %s", prettytime(times[$n]))
uxz_title = @lift @sprintf("u(x, z, t) at y=0 m and t = %s", prettytime(times[$n]))

fig = Figure(resolution = (800, 800))

ax_B = Axis(fig[1, 4];
            xlabel = "Buoyancy (m s⁻²)",
            ylabel = "z (m)")

ax_U = Axis(fig[2, 4];
            xlabel = "Velocities (m s⁻¹)",
            ylabel = "z (m)",
            limits = ((-0.07, 0.07), nothing))

ax_fluxes = Axis(fig[3, 4];
                 xlabel = "Momentum fluxes (m² s⁻²)",
                 ylabel = "z (m)",
                 limits = ((-3.5e-5, 3.5e-5), nothing))

ax_wxy = Axis(fig[1, 1:2];
              xlabel = "x (m)",
              ylabel = "y (m)",
              aspect = DataAspect(),
              limits = ((0, grid.Lx), (0, grid.Ly)),
              title = wxy_title)

ax_wxz = Axis(fig[2, 1:2];
              xlabel = "x (m)",
              ylabel = "z (m)",
              aspect = AxisAspect(2),
              limits = ((0, grid.Lx), (-grid.Lz, 0)),
              title = wxz_title)

ax_uxz = Axis(fig[3, 1:2];
              xlabel = "x (m)",
              ylabel = "z (m)",
              aspect = AxisAspect(2),
              limits = ((0, grid.Lx), (-grid.Lz, 0)),
              title = uxz_title)

nothing #hide

# We use Makie's `Observable` to animate the data. To dive into how `Observable`s work we
# refer to [Makie.jl's Documentation](https://makie.juliaplots.org/stable/documentation/nodes/index.html).

wₙ = @lift time_series.w[$n]
uₙ = @lift time_series.u[$n]
Bₙ = @lift time_series.B[$n][1, 1, :]
Uₙ = @lift time_series.U[$n][1, 1, :]
Vₙ = @lift time_series.V[$n][1, 1, :]
wuₙ = @lift time_series.wu[$n][1, 1, :]
wvₙ = @lift time_series.wv[$n][1, 1, :]

k = searchsortedfirst(grid.zᵃᵃᶠ[:], -8)
wxyₙ = @lift collect(time_series.w[$n][:, :, k])
wxzₙ = @lift collect(time_series.w[$n][:, 1, :])
uxzₙ = @lift collect(time_series.u[$n][:, 1, :])

wlims = (-0.03, 0.03)
ulims = (-0.05, 0.05)

lines!(ax_B, Bₙ, zu)

lines!(ax_U, Uₙ, zu; label = L"\bar{u}")
lines!(ax_U, Vₙ, zu; label = L"\bar{v}")
axislegend(ax_U; position = :rb)

lines!(ax_fluxes, wuₙ, zw; label = L"mean $wu$")
lines!(ax_fluxes, wvₙ, zw; label = L"mean $wv$")
axislegend(ax_fluxes; position = :rb)

hm_wxy = heatmap!(ax_wxy, xw, yw, wxyₙ;
                  colorrange = wlims,
                  colormap = :balance)

Colorbar(fig[1, 3], hm_wxy; label = "m s⁻¹")

hm_wxz = heatmap!(ax_wxz, xw, zw, wxzₙ;
                  colorrange = wlims,
                  colormap = :balance)

Colorbar(fig[2, 3], hm_wxz; label = "m s⁻¹")

ax_uxz = heatmap!(ax_uxz, xu, zu, uxzₙ;
                  colorrange = ulims,
                  colormap = :balance)

Colorbar(fig[3, 3], ax_uxz; label = "m s⁻¹")

# And, finally, we record a movie.

frames = 1:length(times)

record(fig, "langmuir_turbulence.mp4", frames, framerate=8) do i
       @info "Plotting frame $i of $(frames[end])..."
       n[] = i
end
nothing #hide

# ![](langmuir_turbulence.mp4)


text:
#
println("Af) Test of the frequency-independent Breit interaction for an internally generated chlorine-like multiplet.")
#
wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=JAC.Radial.Grid(true), nuclearModel=Nuclear.Model(54.),
                        configs=[Configuration("[Ne] 3s^2 3p^5"), Configuration("[Ne] 3s 3p^6"), Configuration("[Ne] 3s^2 3p^4 3d")],
                        asfSettings=AsfSettings(true, false, Basics.DFSField(), "hydrogenic", Dict{Subshell, Orbital}(), [1],    40, 1.0e-6, JAC.Subshell[], JAC.Subshell[],
                                                true, false, NoneQed(), LSjjSettings(true),
                                                false, [1,2,3,4], false, JAC.LevelSymmetry[] )  )

wb = perform(wa)


text:
module ParameterSchedulers

using Base.Iterators
using Flux
using InfiniteArrays: OneToInf

export Sequence, Loop,
       Step, Exp, Poly, Inv,
       Triangle, TriangleDecay2, TriangleExp,
       Sin, SinDecay2, SinExp,
       Cos

include("decay.jl")
include("cyclic.jl")
include("complex.jl")

# TODO
# Remove this once Optimisers.jl has support
# for schedules + optimizers
"""
    Scheduler{T, O, F}(schedule::AbstractSchedule, opt, update_func)
    Scheduler(schedule, opt; update_func = (o, s) -> (o.eta = s))

Wrap a `schedule` and `opt` together with a `Scheduler`.
The `schedule` is iterated on every call to
[`Flux.apply!`](https://github.com/FluxML/Flux.jl/blob/master/src/optimise/optimisers.jl).
The `Scheduler` can be used anywhere a Flux optimizer is used.

By default, the learning rate (i.e. `opt.eta`) is scheduled.
Set `update_func = (opt, schedule_val) -> ...` to schedule an alternate field.
If `opt` does not have a field `eta`, then there is no default behavior
(you must manually set `update_func`).

# Arguments
- `schedule`: the schedule to use
- `opt`: a Flux optimizer
- `update_func`: a mutating function of with inputs `(optim, param)`
                 that mutates `optim`'s fields based on the current `param` value

# Examples
```julia
# cosine annealing schedule for Descent
julia> s = Cos(λ0 = 0.1, λ1 = 0.8, period = 10);

julia> opt = Scheduler(s, Descent())
Scheduler(Cos{Float64,Int64}(0.1, 0.8, 10), Descent(0.1))

# schedule the momentum term of Momentum
julia> opt = Scheduler(s, Momentum(); update_func = (o, s) -> o.rho = s)
Scheduler(Cos{Float64,Int64}(0.1, 0.8, 10), Momentum(0.01, 0.9, IdDict{Any,Any}()))
```
"""
mutable struct Scheduler{T, O, F}
    state::IdDict{Any, Int}
    schedule::T
    optim::O
    update_func::F
end
Scheduler(schedule, opt, update_func) =
    Scheduler(IdDict{Any, Int}(), schedule, opt, update_func)

Base.show(io::IO, s::Scheduler) =
    print(io, "Scheduler(", s.schedule, ", ", s.optim, ")")

function Flux.Optimise.apply!(opt::Scheduler, x, Δ)
    # get iteration
    t = get!(opt.state, x, 1)
    opt.state[x] = t + 1

    # set param
    opt.update_func(opt.optim, opt.schedule(t))

    # do normal apply
    return Flux.Optimise.apply!(opt.optim, x, Δ)
end

for Opt in (Descent, Momentum, Nesterov, RMSProp,
            ADAM, RADAM, AdaMax, OADAM, ADAGrad,
            ADADelta, AMSGrad, NADAM, AdaBelief)
    @eval begin
        Scheduler(schedule, opt::$Opt; update_func = (o, s) -> (o.eta = s)) =
            Scheduler(schedule, opt, update_func)
    end
end

end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

using Test, Tensors

let dtime = 0.25,
    R0 = 100.0,
    parameters = ChabocheParameterState(E=200.0e3,
                                        nu=0.3,
                                        R0=R0,
                                        Kn=100.0,
                                        nn=10.0,
                                        C1=10000.0,
                                        D1=100.0,
                                        C2=50000.0,
                                        D2=1000.0,
                                        Q=0.0,
                                        b=0.1),
    material = Chaboche(parameters=parameters),
    times = [material.drivers.time],
    stresses = [copy(tovoigt(material.variables.stress))],
    strains = [copy(tovoigt(material.drivers.strain; offdiagscale=2.0))],
    cumeqs = [copy(material.variables.cumeq)],
    tostrain(vec) = fromvoigt(Symm2, vec; offdiagscale=2.0),
    tostress(vec) = fromvoigt(Symm2, vec),
    uniaxial_stress(sigma) = tostress([sigma, 0, 0, 0, 0, 0]),
    stresses_expected = [uniaxial_stress(R0 / 2),
                         uniaxial_stress(R0),
                         uniaxial_stress(1.5 * R0),
                         uniaxial_stress(1.5 * R0),
                         uniaxial_stress(R0),
                         uniaxial_stress(-R0)],
    dstress = R0 / 2,
    dstresses11 = dstress*[1.0, 1.0, 1.0, 0.0, -1.0, -4.0]
    dtimes = [dtime, dtime, dtime, 1e3, dtime, 1e3]

    for i in 1:length(dtimes)
        dstress11 = dstresses11[i]
        dtime = dtimes[i]
        stress_driven_uniaxial_increment!(material, dstress11, dtime)
        update_material!(material)
        push!(times, material.drivers.time)
        push!(stresses, copy(tovoigt(material.variables.stress)))
        push!(strains, copy(tovoigt(material.drivers.strain; offdiagscale=2.0)))
        push!(cumeqs, copy(material.variables.cumeq))
        @test isapprox(material.variables.stress, stresses_expected[i]; atol=1e-4)
    end

    # Plastic creep should have occurred at the portion of the test
    # where the stress was held at 1.5*R0.
    dstrain_creep = strains[5] - strains[4]
    # von Mises material, so the plastic nu = 0.5.
    @test isapprox(dstrain_creep[2], -dstrain_creep[1]*0.5; atol=1e-4)  # ε22 = -0.5 ε11
    @test isapprox(dstrain_creep[3], -dstrain_creep[1]*0.5; atol=1e-4)  # ε33 = -0.5 ε11

    dcumeq = cumeqs[end] - cumeqs[end-1]
    @test dcumeq > 0
end


text:
## Settings Struct

"Store the settings required for performing EEMD."
struct EEMDSetting
    "EMD Setting"
    emd_setting::EMDSetting
    "Ensemble Size"
    ensemble_size::Int64
    "Noise Strength"
    noise_strength::Float64
    "Seed"
    rng_seed::Int

    function EEMDSetting(n::Int64, num_siftings::Int64, s_num::Int64, m::Int64,
                       ensemble_size::Int64, noise_strength::Float64, rng_seed::Int=0)
	    if ensemble_size < 1
            throw(DomainError("Invalid Ensemble Size for EEMD"))
	    elseif noise_strength < 0
            throw(DomainError("Invalid Noise Strength for EMD"))
	    elseif ensemble_size == 1 && noise_strength > 0
            throw(DomainError("Noise Added to EMD"))
    	elseif ensemble_size > 1 && noise_strength == 0
            throw(DomainError("No Noise Added to EEMD"))
        end

        new(EMDSetting(n, num_siftings, s_num, m), ensemble_size, noise_strength, rng_seed)
    end
end

## Main Function

"""
    eemd(input::Vector{Float64}, s::EEMDSetting)

Return the IMFs computed by EEMD given the settings.
"""
function eemd(input::Vector{Float64}, s::EEMDSetting)
    n = length(input)

    if (n == 0)
        throw(DomainError("Invalid size of input"))
    end

    if s.noise_strength != 0
        noise_sigma = std(input)*s.noise_strength
    else
        noise_sigma = zero(s.noise_strength)
    end

    output = zeros(n, s.emd_setting.m)
    rng = Normal(0.0, 1.0)

    for en_i=1:s.ensemble_size
        inp = zeros(n)

        if s.noise_strength == 0.0
            inp = deepcopy(input)
	    else
            Random.seed!(s.rng_seed + en_i)
            rng = Normal(0.0, noise_sigma)

            for i=1:n
                inp[i] += input[i] + rand(rng)
		    end
	    end

        output += emd(inp, s.emd_setting)
    end

    if s.ensemble_size != 1
        output *= 1.0/float(s.ensemble_size)
    end

    return output
end


text:
function inpolygon{T<:Number}(x:: T, y:: T, vx:: Vector{T}, vy:: Vector{T})
    @assert length(vx) == length(vy)
    c = false
    j = length(vx)
    @inbounds for i=1:length(vx)
        if (((vy[i] <= y && y < vy[j]) ||
            (vy[j] <= y && y < vy[i])) &&
            (x < (vx[j] - vx[i]) * (y - vy[i]) / (vy[j] - vy[i]) + vx[i]))
            c = !c
        end
        j = i
    end
    return c
end

function fillpoly!{T,P<:Number}(M::Matrix{T}, px::Vector{P}, py::Vector{P}, value::T)
    @assert length(px) == length(py)
    left, right = round(Int, minimum(px)), round(Int, maximum(px))
    top, bottom = round(Int, minimum(py)), round(Int, maximum(py))
    for x=left:right
        ys = Set{Int64}()
        j = length(px)
        for i=1:length(px)
            if (px[i] <= x <= px[j]) || (px[j] <= x <= px[i])
                # special case: adding the whole cut to ys
                if abs(px[i] - px[j]) <= eps(P)
                    push!(ys, round(Int, py[i]))
                    push!(ys, round(Int, py[j]))
                else
                    y = py[i] + (x - px[i]) / (px[j] - px[i]) * (py[j] - py[i])
                    push!(ys, round(Int, y))
                end
            end
            j = i
        end
        ys = sort([y for y in ys])
        # if there's an odd number of intersection points, add one imeginary point
        if length(ys) % 2 == 1
            push!(ys, ys[end])
        end
        for i=1:2:length(ys)
            M[ys[i]:ys[i+1], x] = value  # <-- bounds error here!
        end
    end
    return M
end

function poly2mask{P}(px::Vector{P}, py::Vector{P}, m::Int, n::Int)
    mask = zeros(Bool, m, n)
    fillpoly!(mask, px, py, true)
end


text:
pmid = 11748933
pmid_list = [24008025, 24170597]

@testset "Export EndNote Citations" begin

    println("-----------------------------------------")
    println("       Export EndNote Citations")

    citation_type="endnote"
    output_file1="11748933.enw"
    output_file2 = "pmid_list.enw"

    export_citation(pmid_list, citation_type, output_file2)
    export_citation(pmid, citation_type, output_file1)

    #Read back the files and do minimal testing
    ref_lines=[]
    open("./11748933.enw") do f
       ref_lines = readlines(f)
    end

    @test ref_lines[1]== "%0 Journal Article"

    ref_lines=[]
    open("./pmid_list.enw") do f
       ref_lines = readlines(f)
    end

    @test ref_lines[1]== "%0 Journal Article"
    @test ref_lines[25]== "%0 Journal Article"

    # remove temp files
    if isfile(output_file1)
        rm(output_file1)
    end
    if isfile(output_file2)
        rm(output_file2)
    end
end

@testset "Export BibTEX Citation" begin
    println("-----------------------------------------")
    println("       Export BibTEX Citations")
    citation_type="bibtex"
    output_file1="11748933.bibtex"
    output_file2 = "pmid_list.bibtex"

    export_citation(pmid, citation_type, output_file1)
    export_citation(pmid_list, citation_type, output_file2)

    #Read back the files and do minimal testing
    ref_lines=[]
    open("./11748933.bibtex") do f
       ref_lines = readlines(f)
    end

    @test ref_lines[1]== "@article {PMID:11748933,"

    ref_lines=[]
    open("./pmid_list.bibtex") do f
       ref_lines = readlines(f)
    end

    @test ref_lines[1]== "@article {PMID:24008025,"
    @test ref_lines[13]== "@article {PMID:24170597,"

    # remove temp files
    if isfile(output_file1)
        rm(output_file1)
    end
    if isfile(output_file2)
        rm(output_file2)
    end
end


text:
### ============================================================================
### Lasso but on linear programming and Quantile Regression
### ============================================================================


text:
function EventFilter(arg0::ODEEventHandler, arg1::FilterType)
    return EventFilter((ODEEventHandler, FilterType), arg0, arg1)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function event_occurred(obj::EventFilter, arg0::ODEStateAndDerivative, arg1::jboolean)
    return jcall(obj, "eventOccurred", Action, (ODEStateAndDerivative, jboolean), arg0, arg1)
end

function g(obj::EventFilter, arg0::ODEStateAndDerivative)
    return jcall(obj, "g", jdouble, (ODEStateAndDerivative,), arg0)
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function init(obj::EventFilter, arg0::ODEStateAndDerivative, arg1::jdouble)
    return jcall(obj, "init", void, (ODEStateAndDerivative, jdouble), arg0, arg1)
end

function reset_state(obj::EventFilter, arg0::ODEStateAndDerivative)
    return jcall(obj, "resetState", ODEState, (ODEStateAndDerivative,), arg0)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end


text:
export IdentitySkip

struct IdentitySkip
   inner
   activation
end

(m::IdentitySkip)(x) = m.activation.(m.inner(x) .+ x)


text:
# Loads trials and computes the weights for the given number of robots

using SubmodularMaximization
using PyPlot
using Base.Iterators
using Base.Threads
using JLD2
using Statistics
using Printf
using Random
using RosDataProcess
using HistogramFilters

close("all")

experiment_name = "weights_by_number_of_robots"
data_folder = "./data"

reprocess = false

# Values taken from the entropy_by_solver script
num_robots = 8:8:40
trials = 1:20
solver_ind = 4 # dist-4

# Drop the first n and keep every few after
trial_steps = 20:100

all_tests = product(num_robots, trials)
all_configurations = num_robots

#
# Load and preprocess prior entropy results
#

solver_results_file = "./data/entropy_by_solver.jld2"

let
  @load solver_results_file results
  global solver_results = results
end

# Pull the relevant trials from the dataset
weights_trials = Dict(map(all_tests) do test
  test => solver_results[solver_ind, test...]
end)

#
# Compute weights
#

function trial_fun(x)
  num_robots, trial = x

  # The nomenclature here is weird, but we will only be using a part of
  # the trial and will implement that at computation time.
  data, configs = weights_trials[x...]
  trial_data = data[trial_steps]

  # Iterate through the trial and compute weights
  trial_weights = map(trial_data) do step_data
    robot_states = step_data.robot_states
    target_filters = step_data.target_filters

    problem = MultiRobotTargetTrackingProblem(robot_states, target_filters,
                                              configs)

    compute_weight_matrix(problem)
  end

  (trial_weights=trial_weights, trial_data=trial_data, configs=configs)
end
function print_summary(x)
  num_robots, trial = x

  println("Num. Robots: ", num_robots,
          " Trial: ", trial)
end

@time results = run_experiments(all_tests,
                                trial_fun=trial_fun,
                                print_summary=print_summary,
                                experiment_name=experiment_name,
                                data_folder=data_folder,
                                reprocess=reprocess
                               )

#
# Plot weights by objective values
#

normalized_weights = map(num_robots) do num_robots
  trial_weights = map(trials) do trial
    trial = results[num_robots, trial]

    trial_weights = trial.trial_weights
    objectives = map(x->x.objective, trial.trial_data)

    total_weights = map(total_weight, trial_weights)

    total_weights ./ objectives
  end

  vcat(trial_weights...)
end

#boxplot([normalized_weights[n] for n in num_robots],
        #notch=false)

weights_series = TimeSeries(num_robots, hcat(normalized_weights...)')

plot_trials(weights_series, mean=true, marker=".", markersize=25, linewidth=3)

title("Target weights")
ylabel("Redundancy over obj.")
xlabel("Number of robots")

save_fig("fig", "weights_by_number_of_robots")

#
# Weights per num robots
#

figure()

weight_per_robot = map(num_robots) do num_robots
  trial_weights = map(trials) do trial
    trial = results[num_robots, trial]

    trial_weights = trial.trial_weights

    total_weights = map(total_weight, trial_weights)

    total_weights / num_robots
  end

  vcat(trial_weights...)
end

boxplot(weight_per_robot, notch=false)

title("Weights per num. robots")
ylabel("Redundancy per robot")
xlabel("Number of robots")
xticks(1:length(num_robots), map(string, num_robots))

save_fig("fig", "weights_per_num_robots")

#
# Plot objective values
#

figure()

objective_values = Dict(map(num_robots) do num_robots
  trial_objectives = map(trials) do trial
    trial = results[num_robots, trial]

    objectives = map(x->x.objective, trial.trial_data)

    objectives / num_robots
  end

  num_robots => vcat(trial_objectives...)
end)

boxplot([objective_values[n] for n in num_robots],
        notch=false)

title("Objective values")
ylabel("Objective value per robot")
xlabel("Number of robots")
xticks(1:length(num_robots), map(string, num_robots))

save_fig("fig", "objective_values_by_number_of_robots")

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
#! /usr/bin/env julia
using PipeLayout
using PipeLayout.GndStr

const ACCOUNT = "gas"

if length(ARGS) != 4
    println("usage: run_all.jl <config> <instlist> <partition> <results>")
    exit(1)
end

CONFIG    = abspath(ARGS[1])
INSTLIST  = ARGS[2]
PARTITION = ARGS[3]
RESULTS   = abspath(ARGS[4])

stem(path) = split(basename(path), ".")[1]

# location of executable
EXEPATH = abspath("run.jl")

# location of instances
INDIR = abspath(dirname(INSTLIST))

# create dir for results
OUTDIR = joinpath(RESULTS, stem(INSTLIST), stem(CONFIG))
mkpath(OUTDIR)

# conservative timelimit (in minutes)
include(CONFIG) # creates solver
@assert solver.timelimit < Inf
timelimit = round(Int, 2 * solver.timelimit / 60)

"submit a job to SLURM"
function submit(key)
    OUT = abspath("$key.log")
    options = ["--account=$ACCOUNT",
               "--partition=$PARTITION",
               "--cpus-per-task=1",
               "--time=$timelimit",
               "--signal=B:INT",
               "--output=$OUT",
               "--error=$OUT"]
    job = [EXEPATH, CONFIG, joinpath(INDIR, key)]
    run(`sbatch $options $job`)
end

# one job per instance, submit from result dir
open(INSTLIST) do f
    for line in eachline(f)
        key = strip(line)
        path = joinpath(OUTDIR, key)

        mkpath(path)
        cd(path)
        submit(key)
    end
end


text:
type SkewTDist <: ContinuousUnivariateDistribution
    ξ::Float64  # Location
    ω::Float64  # Scale
    α::Float64  # Skewness
    df::Float64 # Degrees of freedom
end

skewtdistpdf(α::Real, df::Real, z::Real) = 2.0* _t₁(z, df) * _T₁(α * z * sqrt((df + 1)/(z^2 + df)), df + 1)

function pdf(dist::SkewTDist, x::Real)
    z = (x - dist.ξ)/dist.ω
    skewtdistpdf(dist.α, dist.df, z)/dist.ω
end

function rand(dist::SkewTDist)
    chisqd = Chisq(dist.df)
    w = rand(chisqd)/dist.df
    Ω = Array(Float64, 1,1)
    Ω[1,1] = dist.ω^2
    sndist = MvSkewNormal(Ω, [dist.α])
    x = rand(sndist)[1]
    return dist.ξ + x/sqrt(w)
    ## Ω = Array(Float64, 1,1)
    ## Ω[1,1] = dist.ω^2
    ## return rand(MvSkewTDist([dist.ξ], Ω, [dist.α], dist.df))[1]
end

function cdf(dist::SkewTDist, x::Real)
    quadgk(t->pdf(dist,t), -Inf, x)[1]
end

function quantile(dist::SkewTDist, β::Float64)
    if dist.α < 0
        β = 1 - β
    end
    a = tdistinvcdf(dist.df, β)
    b = sqrt(fdistinvcdf(1.0, dist.df, β))
    qz = fzero(x->cdf(dist, x) - β, a, b)
end

minimum(dist::SkewTDist) = -Inf
maximum(dist::SkewTDist) = Inf

μ(dist::SkewTDist) = δ(dist.α) * sqrt(dist.df/π) * exp(lgamma(0.5*(dist.df-1.0)) - lgamma(0.5*dist.df))
mean(dist::SkewTDist) = dist.ξ + dist.ω * μ(dist)
var(dist::SkewTDist) = (dist.ω^2) * (dist.df/(dist.df - 2.0)) - dist.ω^2 * μ(dist)^2
dof(dist::SkewTDist) = dist.df


text:
__precompile__(true)

module Microeconometrics

using DataFrames
using Dates
using Format: format
using LinearAlgebra
using Optim
using SparseArrays
using SpecialFunctions: loggamma
using StatsFuns
using StatsModels
    const SM = StatsModels
using SuiteSparse

using Base:      @propagate_inbounds
using StatsBase: AbstractWeights, CoefTable, RegressionModel
using StatsBase: AnalyticWeights, FrequencyWeights, ProbabilityWeights, Weights
using StatsBase: mean, sum

import LinearAlgebra: lmul!
import Statistics:    mean
import StatsBase:     varcorrection
import StatsBase:     fit, coef, coefnames, coeftable, confint, stderror, vcov
import StatsBase:     deviance, loglikelihood, nulldeviance, nullloglikelihood
import StatsBase:     adjr2, dof, dof_residual, nobs, r2
import StatsBase:     fitted, predict, residuals, response
import StatsModels:   terms, termvars, schema, apply_schema, has_schema, modelcols

include("./general/types.jl")
include("./data/corrstructures.jl")
include("./data/microdata.jl")
include("./data/micromodel.jl")
include("./data/weights.jl")
include("./data/utils.jl")
include("./general/etable.jl")
include("./general/statsmodels.jl")
include("./general/utils.jl")
include("./inference/covariance_matrix.jl")
include("./inference/model_comparison.jl")
include("./estimation/ols.jl")
include("./estimation/iv.jl")
include("./estimation/logit.jl")
include("./estimation/probit.jl")
include("./estimation/cloglog.jl")
include("./estimation/gompit.jl")
include("./estimation/poisson.jl")
include("./estimation/ivpoisson.jl")
include("./estimation/mullahy.jl")
include("./estimation/ipw.jl")
include("./estimation/abadie.jl")
include("./estimation/frölichmelly.jl")
include("./estimation/tan.jl")

export
    @micromodel, Microdata,
    Homoscedastic, Heteroscedastic, Clustered, CrossCorrelated,
    OLS, Logit, Probit, Cloglog, Gompit, Poisson,
    IV, IVPoisson, Mullahy,
    IPW, Abadie, FrölichMelly, Tan

export
    fit, coef, confint, stderror, vcov,
    deviance, loglikelihood, nulldeviance, nullloglikelihood,
    adjr2, dof, dof_residual, nobs, r2,
    fitted, predict, residuals, response,
    coefnames, coeftable, etable,
    chow_test, hausman_test, pval, tstat

end # module

max_stars_repo_name: lbittarello/Microeconometrics

max_issues_repo_name: lbittarello/Microeconometrics

max_forks_repo_name: lbittarello/Microeconometrics