text:
#   This file is part of Similitude.jl
#   It is licensed under the MIT license
#   Similitude Copyright (C) 2022 Michael Reed
#       _           _                         _
#      | |         | |                       | |
#   ___| |__   __ _| | ___ __ __ ___   ____ _| | __ _
#  / __| '_ \ / _` | |/ / '__/ _` \ \ / / _` | |/ _` |
# | (__| | | | (_| |   <| | | (_| |\ V / (_| | | (_| |
#  \___|_| |_|\__,_|_|\_\_|  \__,_| \_/ \__,_|_|\__,_|
#
#   https://github.com/chakravala
#   https://crucialflow.com

export AbelianGroup, Dimension, 𝟙, normal
export UnitSystems, Quantity, Group, LogGroup, ExpGroup
export universe, Universe, unitname, normal, logdb, expdb, dB

for i ∈ 1:dims
    @eval export $(Symbol(isq[i]))
end

for unit ∈ Dimensionless
    @eval @pure $unit(C::Coupling) = UnitSystems.$unit(C)
    @eval @pure $unit(U::UnitSystem) = UnitSystems.$unit(universe(U))
end

for u ∈ (Constants...,Physics...)
    u∉(:permeability,:gaussgravitation) && @eval const $u = UnitSystems.$u(SI)
end

const hyperfine = SI2019(ΔνCs,inv(T))
const hubble = Hubble(𝟏,inv(T))
const cosmological = 𝟑*ΩΛ*(hubble/lightspeed(Hubble))^2
const eddington = Cosmological(𝟏,M)(QCD)
const solarmass = IAU(𝟏,M)
const earthmass = Metric(GME/G,M)(IAU)
const jupitermass = Metric(GMJ/G,M)(IAU)
const lunarmass = earthmass/μE☾
const gforce = English(𝟏,specificforce)
const atmosphere = Metric(atm,pressure)
const loschmidt = atmosphere(SI2019)/SI2019(T₀,Θ)/boltzmann(SI2019)
const amagat = loschmidt(SI2019)/avogadro(SI2019)
const wienwavelength = planck(SI)*lightspeed(SI)/boltzmann(SI)/Constant(4.965114231744276303)
const wienfrequency = Constant(2.821439372122078893)*boltzmann(SI)/planck(SI)
@pure (::typeof(loschmidt))(U::UnitSystem,P=atmosphere(U),T=SI2019(T₀,Θ)(U)) = U(P,pressure)/U(T,Θ)/boltzmann(U)
@pure mechanicalheat(U::UnitSystem) = molargas(U)*U(normal(calorie(Metric)/molargas(Metric)),Θ*N)

# angle

#const radian = MetricEngineering(𝟏,A)
const spatian = MetricSpatian(𝟏,A)
const steradian = MetricEngineering(𝟏,solidangle)
const degree = MetricDegree(𝟏,A)
const squaredegree = MetricDegree(𝟏,solidangle)
const gradian = MetricGradian(𝟏,A)
const bradian = MetricEngineering(τ/𝟐^8,A)
const arcminute = MetricArcminute(𝟏,A)
const arcsecond = MetricArcsecond(𝟏,A)

# length

const meter = Metric(𝟏,L)
const angstrom = hecto*pico*meter
const inch = IPS(𝟏,L)
#const rackunit = foot*𝟕/𝟐^4/𝟑
const foot = English(𝟏,L)
const surveyfoot = Survey(𝟏,L)
const yard = 𝟑*foot
const mile = English(𝟐^5*𝟑*𝟓*𝟏𝟏,L)
const statutemile = Survey(𝟐^5*𝟑*𝟓*𝟏𝟏,L)
const earthradius = sqrt(earthmass(Metric)*gravitation(Metric)/gforce(Metric))
const greatcircle = τ*earthradius
const earthmeter = Meridian(𝟏,L)
const nauticalmile = Nautical(𝟏,L)
const admiraltymile = English(𝟐^6*𝟓*𝟏𝟗,L)
const meridianmile = Metric(𝟐^4*𝟓^5/𝟑^3,L)
const astronomicalunit = IAU(𝟏,L)
const lunardistance = IAUE(𝟏,L)(IAU)
const jupiterdistance = IAUJ(𝟏,L)(IAU)
const parsec = astronomicalunit*(𝟐^7*𝟑^4*𝟓^3/τ)

#time

const second = Metric(𝟏,T)
const minute = (𝟐^2*𝟑*𝟓)*second
const hour = (𝟐^2*𝟑*𝟓)*minute
const day = IAU(𝟏,T)
const year = IAU(aⱼ,T)
const lightyear = year*lightspeed(IAU)
const radarmile = 𝟐*nauticalmile(Metric)/lightspeed(Metric)
const gaussgravitation = sqrt(normal(gravitation(IAU)))*radian(IAU)/day(IAU)
const gaussianyear = turn(IAU)/gaussgravitation
const siderealyear = gaussianyear/√(solarmass+earthmass+lunarmass).v
const gaussianmonth = τ/sqrt(normal(gravitation(IAUE)))*day
const siderealmonth = gaussianmonth/normal(sqrt(earthmass(IAUE)+lunarmass(IAUE)))
const synodicmonth = inv(inv(siderealmonth(IAU))-inv(siderealyear(IAU)))
const jovianyear = τ*sqrt(normal(jupiterdistance^3/solarmass/gravitation(IAU)))*day/normal(sqrt(solarmass+jupitermass))

# area

const barn = Metric((𝟐*𝟓)^-28,area)
const hectare = Metric(hecto*hecto,area)
const acre = MPH(𝟐^-7/𝟓,area)(English)
const surveyacre = Survey(𝟐^3*𝟑^2*𝟓*𝟏𝟏^2,area)
#const township = 𝟐^9*𝟑^2*𝟓*surveyacre
#const footballfield = English(𝟐^8*𝟑^2*𝟓^2,area)

# volume

const gallon = IPS(𝟑*𝟕*𝟏𝟏,volume)
const liter = Metric(milli,volume)
const quart = gallon/𝟐^2
const pint = quart/𝟐
const cup = pint/𝟐
const fluidounce = cup/𝟐^3
const teaspoon = 𝟓*milli*liter
const tablespoon = 𝟑*teaspoon
#const oilbarrel = 𝟐*𝟑*𝟕*gallon

# speed

const bubnoff = meter(Metric)/year(Metric)
const ips = IPS(𝟏,speed)
const fps = British(𝟏,speed)
const fpm = foot(British)/minute(British)
const ms = Metric(𝟏,speed)
const kmh = kilo*meter/hour
const mph = MPH(𝟏,speed)
const knot = Nautical(𝟏,speed)
const mps = mile(MPH)/second(MPH)

# mass

const gram = Metric(milli,M)
const earthgram = Meridian(milli,M)
const kilogram = Metric(𝟏,M)
const tonne = Metric(kilo,M)
const ton = English(𝟐*kilo,M)
const pound = English(𝟏,M)
const ounce = English(𝟐^-4,M)
const grain = milli*pound/𝟕
const slug = British(𝟏,M)
const slinch = IPS(𝟏,M)
const hyl = GravitationalMetric(𝟏,M)

# force

const dyne = Gauss(𝟏,force)
const newton = Metric(𝟏,force)
const poundal = FPS(𝟏,force)
const kilopond = MetricEngineering(𝟏,force)
const poundforce = English(𝟏,F)

# pressure

const psi = IPS(𝟏,pressure)
const bar = Metric(hecto*kilo,pressure)
const barye = Gauss(𝟏,pressure)
const pascal = Metric(𝟏,pressure)
const technicalatmosphere = kilopond/(centi*meter(ME))^2
const inchmercury = Metric(inv(inHg),pressure)
const torr = Metric(atm/𝟐^3/𝟓/𝟏𝟗,pressure)

# energy

const erg = Gauss(𝟏,energy)
const joule = Metric(𝟏,energy)
const footpound = poundforce*foot
const meancalorie = InternationalMean(𝟐^2*𝟓*𝟑^2/𝟒𝟑,energy)(Metric)
const kilocalorie = International(𝟐^5*𝟓^4*𝟑^2/𝟒𝟑,energy)(Metric)
const calorie = kilocalorie*milli
const earthcalorie = mechanicalheat(Meridian)
const thermalunit = mechanicalheat(English)
const tontnt = giga*calorie(Metric)
const gasgallon = 𝟐*𝟑*𝟏𝟗*kilo*thermalunit(Metric)

# power

const watt = Metric(𝟏,power)
const horsepower = British(𝟐*𝟓^2*𝟏𝟏,power)
const horsepowerwatt = British(𝟐^4*𝟑^3/𝟓*τ,power)
const horsepowermetric = GM(𝟑*𝟓^2,power)
const tonsrefrigeration = thermalunit(Metric)/Metric(𝟑/𝟐/𝟓,T)
const boilerhorsepower = Constant(1339)/Metric(𝟐^4*𝟑^2,T)*thermalunit(Metric)
const electricalhorsepower = Metric(Constant(746),power)

# electromagnetic

const coulomb = Metric(𝟏,Q)
const ampere = Metric(𝟏,current)
const volt = Metric(𝟏,electricpotential)
const henry = Metric(𝟏,inductance)
const ohm = Metric(𝟏,resistance)
const siemens = Metric(𝟏,conductance)
const farad = Metric(𝟏,capacitance)
const weber = Metric(𝟏,magneticflux)
const tesla = Metric(𝟏,magneticfluxdensity)
const statcoulomb = ESU(𝟏,Q)
const statampere = ESU(𝟏,current)
const statvolt = ESU(𝟏,electricpotential)
const stathenry = ESU(𝟏,inductance)
const statohm = ESU(𝟏,resistance)
const statmho = ESU(𝟏,conductance)
const statfarad = ESU(𝟏,capacitance)
const statweber = ESU(𝟏,magneticflux)
const stattesla = ESU(𝟏,magneticfluxdensity)
const abcoulomb = EMU(𝟏,Q)
const abampere = EMU(𝟏,current)
const abvolt = EMU(𝟏,electricpotential)
const abhenry = EMU(𝟏,inductance)
const abohm = EMU(𝟏,resistance)
const abmho = EMU(𝟏,conductance)
const abfarad = EMU(𝟏,capacitance)
const maxwell = EMU(𝟏,magneticflux)
const gauss = EMU(𝟏,magneticfluxdensity)
const oersted = EMU(𝟏,magneticfield)
const gilbert = EMU(𝟏/𝟐/τ,current/A)
const earthcoulomb = Meridian(𝟏,Q)
const electronvolt = elementarycharge(SI2019)*SI2019(𝟏,electricpotential)

# temperature

#const freezing = Metric(T₀-milli,Θ)
const boiling = Metric(T₀+Constant(99.9839),Θ)
const sealevel = Metric(T₀+𝟑*𝟓,Θ)
const kelvin = Metric(𝟏,Θ)
const celsius = Metric(T₀,Θ)
const rankine = English(𝟏,Θ)
const fahrenheit = English(Constant(459.67),Θ)
#const delisle = Metric(𝟐/𝟑,Θ)
#const reaumur = Metric(𝟓/𝟐^2,Θ)

# mole

const mole = Metric(𝟏,N)
const earthmole = Meridian(𝟏,N)
const poundmole = English(𝟏,N)
const slugmole = British(𝟏,N)
const slinchmole = IPS(𝟏,N)

# photometric

const lumen = Metric(𝟏,luminousflux)
const candela = Metric(𝟏,luminousintensity)
const lux = Metric(𝟏,illuminance)
const phot = Gauss(𝟏,illuminance)
const footcandle = English(𝟏,illuminance)
const nit = Metric(𝟏,luminance)
const apostilb = Metric(𝟐/τ,luminance)
const stilb = Gauss(𝟏,luminance)
const lambert = Gauss(𝟐/τ,luminance)
const footlambert = English(𝟐/τ,luminance)
const bril = centi*nano*lambert
const talbot = Metric(𝟏,luminousenergy)
const lumerg = Gauss(centi^2*milli,luminousenergy)
const rayleigh = Metric(deka*giga,photonirradiance)
const flick = Metric(deka*giga,radiance/L)

@pure neper(U::UnitSystem) = U(𝟏,log(𝟙))
@pure bel(U::UnitSystem) = U(𝟏,log10(𝟙))
@pure decibel(U::UnitSystem) = U(𝟏,dB(𝟙))
const hertz = inv(second(Metric))
const apm = inv(minute(Metric))
const rpm = turn(Metric)/minute(Metric)
#const rpd = turn(Metric)/day(Metric)
const galileo = Gauss(𝟏,specificforce)
const eotvos = Gauss(nano,specificforce/L)
const poise = Gauss(𝟏,viscosity)
const reyn = IPS(𝟏,viscosity)
const diopter = Metric(𝟏,angularwavenumber)
const kayser = Gauss(𝟏,wavenumber)
const darcy = Gauss(milli/atm,area)
const stokes = Gauss(𝟏,diffusivity)
const katal = Metric(𝟏,catalysis)
const mpge = mile(Metric)/gasgallon(Metric)
const curie = Constant(37)*giga*hertz
const sievert = Metric(𝟏,energy/M)
#const rem = centi*sievert
const roentgen = ESU(𝟏,chargedensity)(Metric)/Metric(Constant(1.293),density)
const rayl = Metric(𝟏,specificimpedance)
const langley = calorie(Metric)/(centi*meter(Metric))^2
const jansky = Metric((𝟐*𝟓)^-26,fluence)
const solarflux = hecto*hecto*jansky

for CAL ∈ (:calₜₕ,:cal₄,:cal₁₀,:cal₂₀,:calₘ,:calᵢₜ)
    KCAL = Symbol(:k,CAL)
    @eval const $CAL = SI(UnitSystems.$CAL,energy)
    @eval const $KCAL = SI(UnitSystems.$KCAL,energy)
end

#const H0 = hubble

evaldim(::typeof(angle)) = A
evaldim(::typeof(length)) = L
evaldim(::typeof(time)) = T
evaldim(::typeof(molarmass)) = M/N
evaldim(::typeof(luminousefficacy)) = T*J/F/L
evaldim(::typeof(UnitSystems.solidangle)) = A^2
evaldim(unit::Dimension) = unit
evaldim(unit::Symbol) = evaldim(eval(unit))
evaldim(unit::Symbol,U) = evaldim(evaldim(unit),U)
evaldim(unit,U) = normal(U)(unit)

(::typeof(molarmass))(U::UnitSystem,S::UnitSystem) = (M/N)(U,S)
(::typeof(luminousefficacy))(U::UnitSystem,S::UnitSystem) = (T*J/F/L)(U,S)

dimlist(text,dim) = "$text : [$(evaldim(dim))], [$(evaldim(dim,British))], [$(evaldim(dim,Metric))], [$(evaldim(dim,EMU))], [$(evaldim(dim,ESU))]"
dimlist(U) = join(dimtext(normal(U)) == isq ? isodim.(Ref(U),usq) : unitdim.(Ref(U),usq), ", ")
isodim(U,D) = (UD = U(D); D==UD ? "$D" : "$D=$UD")
unitdim(U,D) = (io = IOBuffer(); showgroup(io,U(D),U); "$D=$(String(take!(io)))")

convertext(unit,fun) = """
```Julia
$(dimlist(unit,unit))
$unit(U::UnitSystem,S::UnitSystem) = $fun
$unit(v::Real,U::UnitSystem,S::UnitSystem) = v/$unit(U,S)
$(evaldim(unit)(Unified))
```
"""

@pure unitsym(x) = :nonstandard
@pure unitsym(::typeof(𝟙)) = :dimensionless
@pure unitsym(::typeof(A)) = :angle
for unit ∈ Convert
    if unit ∉ (:length,:time,:angle,:molarmass,:luminousefficacy)
        @eval @pure unitsym(::typeof($(eval(unit)))) = $(QuoteNode(unit))
    else
        @eval @pure unitsym(::typeof($(evaldim(unit)))) = $(QuoteNode(unit))
    end
end

function unitext(unit,text)
    dim = Dimension(eval(unit))
    sym = unitsym(dim)
    return """
```Julia
$unit(U::UnitSystem) = $text
$(dimlist(sym,dim))
$(eval(unit)(Unified))
```
"""
end

systext(sys,text) = """
```Julia
$sys = $text
$(dimlist(eval(sys)))
```
"""

# 1,2,3,4, 5, 6, 7,  8,9,10,11
#kB,ħ,𝘤,μ₀,mₑ,Mᵤ,Kcd,A,λ,αL,g₀
# F,M,L,T, Q, Θ, N,  J,A,R, C

function (u::typeof(normal(MetricEngineering)))(d::Group)
    Group(Values(d.v[1],d.v[2],d.v[3],d.v[4],d.v[5],d.v[6],d.v[7],d.v[8],d.v[9],0,0))
end
function (u::typeof(normal(GravitationalMetric)))(d::Group)
    Group(Values(d.v[1]+d.v[2],0,d.v[3]-d.v[2],d.v[4]+2(d.v[2]),d.v[5],d.v[6],d.v[7],d.v[8],0,0,0))
end
function (u::typeof(normal(Metric)))(d::Group)
    Group(Values(0,d.v[1]+d.v[2],d.v[1]+d.v[3],d.v[4]-2(d.v[1]),d.v[5],d.v[6],d.v[7],d.v[8],0,0,0))
end
function (u::typeof(normal(MetricDegree)))(d::Group)
    Group(Values(0,d.v[1]+d.v[2],d.v[1]+d.v[3],d.v[4]-2(d.v[1]),d.v[5],d.v[6],d.v[7],d.v[8],d.v[9],0,0))
end
#=function (u::typeof(normal(Astronomical)))(d::Group)
    Group(Values(d.v[1],0,d.v[3]+3d.v[2],d.v[4]-2d.v[2],d.v[5],d.v[6],d.v[7],d.v[8],d.v[9],0,0))
end=#

function (u::typeof(normal(Gauss)))(d::Group{<:Integer})
    Group(Values(0,d.v[1]+d.v[2]+d.v[5]//2,d.v[1]+d.v[3]+(3//2)*d.v[5]+d.v[11],d.v[4]-2(d.v[1])-d.v[5]-d.v[11],0,d.v[6],d.v[7],d.v[8],0,0,0))
end
function (u::typeof(normal(ESU)))(d::Group{<:Integer})
    Group(Values(0,d.v[1]+d.v[2]+d.v[5]//2,d.v[1]+d.v[3]+(3//2)*d.v[5],d.v[4]-2(d.v[1])-d.v[5],0,d.v[6],d.v[7],d.v[8],0,0,0))
end
function (u::typeof(normal(EMU)))(d::Group{<:Integer})
    Group(Values(0,d.v[1]+d.v[2]+d.v[5]//2,d.v[1]+d.v[3]+d.v[5]//2,d.v[4]-2(d.v[1]),0,d.v[6],d.v[7],d.v[8],0,0,0))
end

function (u::typeof(normal(Gauss)))(d::Group)
    Group(Values(0,d.v[1]+d.v[2]+d.v[5]/2,d.v[1]+d.v[3]+(3/2)*d.v[5]+d.v[11],d.v[4]-2(d.v[1])-d.v[5]-d.v[11],0,d.v[6],d.v[7],d.v[8],0,0,0))
end
function (u::typeof(normal(ESU)))(d::Group)
    Group(Values(0,d.v[1]+d.v[2]+d.v[5]/2,d.v[1]+d.v[3]+(3/2)*d.v[5],d.v[4]-2(d.v[1])-d.v[5],0,d.v[6],d.v[7],d.v[8],0,0,0))
end
function (u::typeof(normal(EMU)))(d::Group)
    Group(Values(0,d.v[1]+d.v[2]+d.v[5]/2,d.v[1]+d.v[3]+d.v[5]/2,d.v[4]-2(d.v[1]),0,d.v[6],d.v[7],d.v[8],0,d.v[10],0))
end

function (u::typeof(normal(Stoney)))(d::Group)
    Group(Values(0,d.v[1]+d.v[2]+d.v[6]+d.v[7],0,d.v[3]+d.v[4]-d.v[1],d.v[5],0,0,d.v[8],0,0,0))
end
function (u::typeof(normal(Electronic)))(d::Group)
    Group(Values(0,0,0,d.v[3]+d.v[4]-d.v[1]-d.v[8],d.v[5],0,0,0,0,0,0))
end
function (u::typeof(normal(QCDoriginal)))(d::Group)
    Group(Values(0,d.v[2]+d.v[6]+d.v[7]+2(d.v[1]+d.v[8])-d.v[3]-d.v[4],0,0,d.v[5],0,0,0,0,0,0))
end
function (u::typeof(normal(Planck)))(d::Group)
    Group(Values(0,d.v[2]+d.v[6]+d.v[7]+2(d.v[1]+d.v[8])-d.v[3]-d.v[4],0,0,0,0,0,0,0,0,0))
end
function (u::typeof(normal(PlanckGauss)))(d::Group)
    Group(Values(0,d.v[2]+d.v[6]+d.v[7]+2(d.v[1]+d.v[8])-d.v[3]-d.v[4],0,0,d.v[5],0,0,0,0,0,0))
end
function (u::typeof(normal(Natural)))(d::Group)
    Group(Values(0,0,0,0,0,0,0,0,0,0,0))
end
function (u::typeof(normal(NaturalGauss)))(d::Group)
    Group(Values(0,0,0,0,d.v[5],0,0,0,0,0,0))
end

function (u::typeof(normal(Rydberg)))(d::Group)
    Group(Values(0,d.v[1]+d.v[2]+d.v[7],d.v[1]+d.v[3],d.v[4]-d.v[6]+2(d.v[8]-d.v[1]),d.v[5],0,0,0,0,0,0))
end
function (u::typeof(normal(Hartree)))(d::Group)
    Group(Values(0,0,d.v[3]+2(d.v[4]-d.v[6])-3(d.v[1])-4(d.v[8]),0,d.v[5],0,0,0,0,0,0))
end
function (u::typeof(normal(Hubble)))(d::Group)
    Group(Values(0,0,0,d.v[3]+d.v[4]-d.v[1]-d.v[8],d.v[5],0,0,0,0,0,0))
end
function (u::typeof(normal(Cosmological)))(d::Group)
    Group(Values(0,d.v[1]+d.v[2]+d.v[6]+d.v[7],0,d.v[3]+d.v[4]-d.v[1],d.v[5],0,0,d.v[8],0,0,0))
end
function (u::typeof(normal(CosmologicalQuantum)))(d::Group)
    Group(Values(0,d.v[2]+d.v[6]+d.v[7]+2(d.v[1]+d.v[8])-d.v[3]-d.v[4],0,0,d.v[5],0,0,0,0,0,0))
end

export @unitdim, @unitgroup

"""
    @unitgroup(U::UnitSystem,S::UnitSystem) -> (u::typeof(normal(U)))(d::Group) = normal(S)(d)

Implements `Group` homomorphism for `U` in terms of existing specification from `S`.
"""
macro unitgroup(U,S)
    :((u::typeof(normal($U)))(d::Group) = normal($S)(d))
end

#=@unitgroup SI2019 Metric
@unitgroup SI1976 Metric
@unitgroup Conventional Metric
@unitgroup CODATA Metric
@unitgroup International Metric
@unitgroup InternationalMean Metric
@unitgroup MTS Metric
@unitgroup KKH Metric
@unitgroup MPH Metric
@unitgroup Nautical Metric
@unitgroup Meridian Metric
@unitgroup FFF Metric
@unitgroup IAU Metric
@unitgroup IAUE Metric
@unitgroup IAUJ Metric=#

@unitgroup MetricTurn MetricDegree
@unitgroup MetricSpatian MetricDegree
@unitgroup MetricGradian MetricDegree
@unitgroup MetricArcminute MetricDegree
@unitgroup MetricArcsecond MetricDegree
@unitgroup LorentzHeaviside Gauss
#@unitgroup Thomson EMU
#@unitgroup Kennelly EMU
@unitgroup Schrodinger Rydberg
@unitgroup QCD Planck
@unitgroup QCDGauss PlanckGauss
#@unitgroup Cosmological Hubble

#@unitgroup SI2019Engineering MetricEngineering
#@unitgroup MeridianEngineering MetricEngineering
#@unitgroup GravitationalSI2019 GravitationalMetric
#@unitgroup GravitationalMeridian GravitationalMetric
@unitgroup British GravitationalMetric
@unitgroup English MetricEngineering
@unitgroup Survey MetricEngineering
@unitgroup IPS GravitationalMetric

"""
    @unitdim(U::UnitSystem,F,M,L,T,Q,Θ,N,J="lm",A="rad")

Specify the `print` output for each base `Dimension` of `U::UnitSystem` with `String` input arguments `force`, `mass`, `length`, `time`, `charge`, `temperature`, `molaramount`, `luminousflux`, `angle`.
```Julia
@unitdim Gauss "gf" "g" "cm" "s" "C" "K" "mol"
@unitdim Metric "kgf" "kg" "m" "s" "C" "K" "mol"
@unitdim British "lb" "slug" "ft" "s" "C" "°R" "slug-mol"
@unitdim IPS "lb" "slinch" "in" "s" "C" "°R" "slinch-mol"
@unitdim FPS "pdl" "lb" "ft" "s" "C" "°R" "lb-mol"
@unitdim English "lbf" "lbm" "ft" "s" "C" "°R" "lb-mol"
@unitdim IAU☉ "M☉f" "M☉" "au" "D" "C" "K" "mol"
```
These standard examples are some of the built-in defaults.
"""
macro unitdim(U,F,M,L,T,Q,Θ,N,J="lm",A="rad",R="",C="")
    :(dimtext(::typeof(normal($U))) = Values($F,$M,$L,$T,$Q,$Θ,$N,$J,$A,$R,$C))
end

@unitdim Metric "kgf" "kg" "m" "s" "C" "K" "mol"
@unitdim Meridian "kegf" "keg" "em" "s" "eC" "K" "eg-mol"
@unitdim MetricTurn "kgf" "kg" "m" "s" "C" "K" "mol" "lm" "τ"
@unitdim MetricSpatian "kgf" "kg" "m" "s" "C" "K" "mol" "lm" "ς"
@unitdim MetricGradian "kgf" "kg" "m" "s" "C" "K" "mol" "lm" "gon"
@unitdim MetricDegree "kgf" "kg" "m" "s" "C" "K" "mol" "lm" "deg"
@unitdim MetricArcminute "kgf" "kg" "m" "s" "C" "K" "mol" "lm" "amin"
@unitdim MetricArcsecond "kgf" "kg" "m" "s" "C" "K" "mol" "lm" "asec"
@unitdim British "lb" "slug" "ft" "s" "C" "°R" "slug-mol"
@unitdim English "lbf" "lbm" "ft" "s" "C" "°R" "lb-mol"
@unitdim IPS "lb" "slinch" "in" "s" "C" "°R" "slinch-mol"
@unitdim FPS "pdl" "lb" "ft" "s" "C" "°R" "lb-mol"
@unitdim Gauss "gf" "g" "cm" "s" "_" "K" "mol"
@unitdim IAU☉ "M☉f" "M☉" "au" "D" "C" "K" "mol"
@unitdim IAUE "MEf" "ME" "LD" "D" "C" "K" "mol"
@unitdim IAUJ "MJf" "MJ" "JD" "D" "C" "K" "mol"
@unitdim MTS "tf" "t" "m" "s" "C" "K" "mol"
@unitdim KKH "kgf" "kg" "km" "h" "C" "K" "mol"
@unitdim MPH "lbf" "lb" "mi" "h" "C" "°R" "lb-mol"
@unitdim Nautical "kegf" "keg" "nm" "h" "eC" "K" "eg-mol"
@unitdim FFF "firf" "fir" "fur" "ftn" "Inf" "°R" "fir-mol"
@unitdim Hartree "F" "M" "a₀" "T" "𝘦" "Θ" "N" "J"
@unitdim QCDoriginal "F" "mₚ" "L" "T" "𝘦" "Θ" "N" "J"
@unitdim QCD "F" "mₚ" "L" "T" "Q" "Θ" "N" "J"
@unitdim QCDGauss "F" "mₚ" "L" "T" "𝘦ₙ" "Θ" "N" "J"
@unitdim PlanckGauss "F" "mP" "L" "T" "𝘦ₙ" "Θ" "N" "J"
@unitdim NaturalGauss "F" "T" "L" "T" "𝘦ₙ" "Θ" "N" "J"

"""
    @unitdim(U::UnitSystem,S::UnitSystem) -> dimtext(::typeof(normal(U))) = dimtext(normal(S))

Specify the `print` output for each base `Dimension` of `U` upon prior existing `S` data.
```Julia
@unitdim EMU Gauss
@unitdim ESU Gauss
@unitdim LorentzHeaviside Gauss
@unitdim SI2019 Metric
@unitdim SI1976 Metric
@unitdim CODATA Metric
@unitdim Conventional Metric
@unitdim International Metric
@unitdim InternationalMean Metric
@unitdim Survey English
```
These standard examples are some of the built-in defaults.
"""
macro unitdim(U,S)
    :(dimtext(::typeof(normal($U))) = dimtext(normal($S)))
end

@unitdim SI2019 Metric
@unitdim SI1976 Metric
@unitdim CODATA Metric
@unitdim Conventional Metric
@unitdim International Metric
@unitdim InternationalMean Metric
@unitdim MetricEngineering Metric
@unitdim GravitationalMetric Metric
#@unitdim SI2019Engineering MetricEngineering
#@unitdim GravitationalSI2019 GravitationalMetric
#@unitdim MeridianEngineering Meridian
#@unitdim GravitationalMeridian Meridian
@unitdim Survey English
@unitdim EMU Gauss
@unitdim ESU Gauss
@unitdim LorentzHeaviside Gauss
#@unitdim Thomson Gauss
#@unitdim Kennelly Metric

"""
    @unitdim(D,U,S) -> showgroup(io::IO,::typeof(U(D)),::typeof(normal(U))) = print(io,S)

Specify the `print` output `S::String` for derived `D::Dimension` in `U::UnitSystem`.
```Julia
@unitdim magneticflux Gauss "Mx"
@unitdim magneticfluxdensity Gauss "G"
@unitdim magneticfield Gauss "Oe"
@unitdim frequency Metric "Hz"
@unitdim force Metric "N"
@unitdim pressure Metric "Pa"
@unitdim energy Metric "J"
@unitdim power Metric "W"
@unitdim mass British "slug"
@unitdim force FPS "pdl"
```
These standard examples are some of the built-in defaults.
"""
macro unitdim(D, U, S)
    :(showgroup(io::IO,::typeof($U($D)),::typeof(normal($U))) = print(io,$S))
end

for U ∈ (:MetricEngineering,:GravitationalMetric)
    @eval begin
        @unitdim frequency $U "Hz"
        @unitdim frequencydrift $U "Hz⋅s⁻¹"
        @unitdim photonirradiance $U "Hz⋅m⁻²"
        @unitdim illuminance $U "lx"
        @unitdim luminousexposure $U "lx⋅s"
        showgroup(io::IO,::typeof(luminance),::typeof(normal($U))) = print(io,"nt")
    end
end
for U ∈ (:MetricEngineering,:English,:Survey)
    @eval @unitdim specificforce $U "g₀"
end
for U ∈ (:Metric, :SI2019, :CODATA, :Conventional, :International, :InternationalMean, :MetricTurn, :MetricSpatian, :MetricGradian, :MetricDegree, :MetricArcminute, :MetricArcsecond)
    @eval begin
        @unitdim frequency $U "Hz"
        @unitdim frequencydrift $U "Hz⋅s⁻¹"
        @unitdim photonirradiance $U "Hz⋅m⁻²"
        @unitdim force $U "N"
        @unitdim inv(force) $U "N⁻¹"
        @unitdim pressure $U "Pa"
        @unitdim compressibility $U "Pa⁻¹"
        @unitdim energy $U "J"
        @unitdim inv(energy) $U "J⁻¹"
        @unitdim power $U "W"
        @unitdim inv(power) $U "W⁻¹"

        @unitdim electricpotential $U "V"
        @unitdim inv(electricpotential) $U "V⁻¹"
        @unitdim capacitance $U "F"
        @unitdim inv(capacitance) $U "F⁻¹"
        @unitdim resistance $U "Ω"
        @unitdim conductance $U "S"
        @unitdim magneticflux $U "Wb"
        @unitdim inv(magneticflux) $U "Hz⋅V⁻¹"
        @unitdim magneticfluxdensity $U "T"
        @unitdim inv(magneticfluxdensity) $U "T⁻¹"
        @unitdim permeance $U "H"
        @unitdim reluctance $U "H⁻¹"

        @unitdim catalysis $U "kat"
        @unitdim molarenergy $U "J⋅mol⁻¹"
        @unitdim molarentropy $U "J⋅K⁻¹mol⁻¹"

        @unitdim luminousflux/power $U "lm⋅W⁻¹"
        @unitdim power/luminousflux $U "W⋅lm⁻¹"
        @unitdim illuminance $U "lx"
        @unitdim luminousexposure $U "lx⋅s"

        @unitdim action*speed $U "J⋅m"
        @unitdim impulse $U "N⋅s"
        @unitdim yank $U "N⋅s⁻¹"
        @unitdim fluence $U "N⋅m⁻¹"
        @unitdim compliance $U "m⋅N⁻¹"

        @unitdim viscosity $U "Pa⋅s"
        @unitdim irradiance $U "W⋅m⁻²"
        @unitdim inv(irradiance) $U "W⁻¹m²"
        @unitdim powerdensity $U "W⋅m⁻³"
        @unitdim spectralexposure $U "J⋅m⁻²⋅Hz⁻¹"
        @unitdim irradiance/Θ^4 $U "W⋅m⁻²K⁻⁴"
        @unitdim pressure/Θ^4 $U "J⋅m⁻³K⁻⁴"
        @unitdim 𝟙/T/Θ $U "Hz⋅K⁻¹"
        @unitdim entropy/Q $U "V⋅K⁻¹"
        @unitdim entropy $U "J⋅K⁻¹"
        @unitdim specificentropy $U "J⋅K⁻¹kg⁻¹"
        @unitdim specificenergy $U "J⋅kg⁻¹"
        @unitdim thermalconductivity $U "W⋅m⁻¹K⁻¹"
        @unitdim thermalconductance $U "W⋅K⁻¹"
        @unitdim thermalresistance $U "K⋅W⁻¹"
        @unitdim thermalresistivity $U "K⋅m⋅W⁻¹"
        @unitdim molarconductivity $U "S⋅m²mol⁻¹"

        @unitdim electricpotential/M $U "V⋅kg⁻¹"
        @unitdim electricflux $U "V⋅m"
        @unitdim electricfield $U "V⋅m⁻¹"
        @unitdim permittivity $U "F⋅m⁻¹"
        @unitdim inv(permittivity) $U "m⋅F⁻¹"
        @unitdim permeability $U "H⋅m⁻¹"
        @unitdim inv(permeability) $U "m⋅H⁻¹"
        @unitdim resistivity $U "Ω⋅m"
        @unitdim conductivity $U "S⋅m⁻¹"
        @unitdim vectorpotential $U "Wb⋅m⁻¹"
        @unitdim magneticmoment $U "Wb⋅m"
        @unitdim mobility $U "m²s⁻¹V⁻¹"
    end
end
for U ∈ (:Metric, :SI2019, :CODATA, :Conventional, :International, :InternationalMean)
    @eval begin
        @unitdim luminousintensity $U "cd"
        showgroup(io::IO,::typeof(luminance),::typeof(normal($U))) = print(io,"nt")
        @unitdim angularmomentum $U "J⋅s"
        @unitdim magneticdipolemoment $U "J⋅T⁻¹"
    end
end
for U ∈ (:MetricTurn,:MetricSpatian,:MetricDegree,:MetricGradian,:MetricArcminute,:MetricArcsecond)
    let u = dimtext(normal(eval(U)))[9]
        @eval begin
            @unitdim angularmomentum $U $("J⋅s⋅$(u)⁻¹")
            @unitdim magneticdipolemoment $U $("J⋅T⁻¹⋅$(u)⁻¹")
            @unitdim photonintensity $U $("Hz⋅$(u)⁻²")
            @unitdim photonradiance $U $("Hz⋅m⁻²⋅$(u)⁻²")
            @unitdim radiance $U $("W⋅m⁻²⋅$(u)⁻²")
            @unitdim radiance*T $U $("W⋅m⁻²⋅$(u)⁻²⋅Hz⁻¹")
            @unitdim radiance/L $U $("W⋅m⁻³⋅$(u)⁻²")
            @unitdim radiantintensity $U "W⋅$(u)⁻²"
            @unitdim radiantintensity*T $U $("W⋅$(u)⁻²⋅Hz⁻¹")
            @unitdim radiantintensity/L $U $("W⋅$(u)⁻²⋅m⁻¹")
        end
    end
end

@unitdim frequency  Meridian "Hz"
@unitdim frequencydrift Meridian "Hz⋅s⁻¹"
@unitdim photonirradiance Meridian "Hz⋅m⁻²"
@unitdim force Meridian "eN"
@unitdim inv(force) Meridian "eN⁻¹"
@unitdim pressure Meridian "ePa"
@unitdim compressibility Meridian "ePa⁻¹"
@unitdim energy Meridian "eJ"
@unitdim inv(energy) Meridian "eJ⁻¹"
@unitdim power Meridian "eW"
@unitdim inv(power) Meridian "eW⁻¹"

@unitdim electricpotential Meridian "eV"
@unitdim inv(electricpotential) Meridian "eV⁻¹"
@unitdim capacitance Meridian "eF"
@unitdim inv(capacitance) Meridian "eF⁻¹"
@unitdim resistance Meridian "eΩ"
@unitdim conductance Meridian "eS"
@unitdim magneticflux Meridian "eWb"
@unitdim inv(magneticflux) Meridian "Hz⋅eV⁻¹"
@unitdim magneticfluxdensity Meridian "eT"
@unitdim inv(magneticfluxdensity) Meridian "eT⁻¹"
@unitdim permeance Meridian "eH"
@unitdim reluctance Meridian "eH⁻¹"

@unitdim catalysis Meridian "ekat"
@unitdim molarenergy Meridian "eJ⋅eg-mol⁻¹"
@unitdim molarentropy Meridian "eJ⋅K⁻¹eg-mol⁻¹"

@unitdim luminousflux/power Meridian "lm⋅eW⁻¹"
@unitdim luminousintensity Meridian "cd"
@unitdim illuminance Meridian "elx"
@unitdim luminousexposure Meridian "lx⋅s"
showgroup(io::IO,::typeof(luminance),::typeof(normal(Meridian))) = print(io,"ent")

@unitdim impulse Meridian "eN⋅s"
@unitdim angularmomentum Meridian "eJ⋅s"
@unitdim action*speed Meridian "eJ⋅em"
@unitdim yank Meridian "eN⋅s⁻¹"
@unitdim fluence Meridian "eN⋅em⁻¹"
@unitdim compliance Meridian "em⋅eN⁻¹"

@unitdim viscosity Meridian "ePa⋅s"
@unitdim irradiance Meridian "eW⋅em⁻²"
@unitdim inv(irradiance) Meridian "eW⁻¹em²"
@unitdim powerdensity Meridian "eW⋅m⁻³"
@unitdim spectralexposure Meridian "eJ⋅em⁻²⋅Hz⁻¹"
@unitdim irradiance/Θ^4 Meridian "eW⋅em⁻²K⁻⁴"
@unitdim pressure/Θ^4 Meridian "eJ⋅em⁻³K⁻⁴"
@unitdim 𝟙/T/Θ Meridian "Hz⋅K⁻¹"
@unitdim entropy/Q Meridian "eV⋅K⁻¹"
@unitdim entropy Meridian "eJ⋅K⁻¹"
@unitdim specificentropy Meridian "eJ⋅K⁻¹keg⁻¹"
@unitdim specificenergy Meridian "eJ⋅keg⁻¹"
@unitdim thermalconductivity Meridian "eW⋅em⁻¹K⁻¹"
@unitdim thermalresistance Meridian "K⋅eW⁻¹"
@unitdim thermalresistivity Meridian "K⋅em⋅eW⁻¹"
@unitdim molarconductivity Meridian "eS⋅em²eg-mol⁻¹"

@unitdim electricpotential/M Meridian "eV⋅kg⁻¹"
@unitdim action*speed/Q Meridian "eV⋅em"
@unitdim electricfield Meridian "eV⋅em⁻¹"
@unitdim permittivity Meridian "eF⋅em⁻¹"
@unitdim inv(permittivity) Meridian "em⋅eF⁻¹"
@unitdim permeability Meridian "eH⋅em⁻¹"
@unitdim inv(permeability) Meridian "em⋅eH⁻¹"
@unitdim resistivity Meridian "eΩ⋅em"
@unitdim conductivity Meridian "eS⋅em⁻¹"
@unitdim magneticdipolemoment Meridian "eJ⋅eT⁻¹"
@unitdim vectorpotential Meridian "eWb⋅em⁻¹"
@unitdim magneticmoment Meridian "eWb⋅em"
@unitdim mobility Meridian "em²s⁻¹eV⁻¹"

for U ∈ (:Gauss, :EMU, :ESU, :LorentzHeaviside)
    @eval begin
        @unitdim volume $U "mL"
        @unitdim numberdensity $U "mL⁻¹"
        @unitdim frequency $U "Hz"
        @unitdim photonirradiance $U "Hz⋅m⁻²"
        @unitdim force $U "dyn"
        @unitdim inv(force) $U "dyn⁻¹"
        @unitdim specificforce $U "gal"
        @unitdim specificforce/L $U "gal⋅cm⁻¹"
        @unitdim pressure $U "Ba"
        @unitdim compressibility $U "Ba⁻¹"
        @unitdim energy $U "erg"
        @unitdim inv(energy) $U "erg⁻¹"
        @unitdim power $U "erg⋅s⁻¹"
        @unitdim inv(power) $U "s⋅erg⁻¹"

        @unitdim catalysis $U "kat"
        @unitdim molarenergy $U "erg⋅mol⁻¹"
        @unitdim molarentropy $U "erg⋅K⁻¹mol⁻¹"

        @unitdim luminousflux/power $U "lm⋅s⋅erg⁻¹"
        @unitdim power/luminousflux $U "erg⋅s⁻¹lm⁻¹"
        @unitdim luminousintensity $U "cd"
        @unitdim illuminance $U "ph"
        showgroup(io::IO,::typeof(luminance),::typeof(normal($U))) = print(io,"sb")

        @unitdim angularmomentum $U "erg⋅s"
        @unitdim action*speed $U "erg⋅cm"
        @unitdim fluence $U "dyn⋅cm⁻¹"
        @unitdim compliance $U "cm⋅dyn⁻¹"
        @unitdim impulse $U "dyn⋅s"
        @unitdim yank $U "dyn⋅s⁻¹"

        @unitdim viscosity $U "P"
        @unitdim diffusivity $U "St"
        @unitdim irradiance $U "erg⋅s⁻¹cm⁻²"
        @unitdim inv(irradiance) $U "erg⁻¹s⋅cm²"
        @unitdim powerdensity $U "erg⋅s⁻¹mL⁻¹"
        @unitdim spectralexposure $U "erg⋅cm⁻²⋅Hz⁻¹"
        @unitdim irradiance/Θ^4 $U "erg⋅s⁻¹cm⁻²K⁻⁴"
        @unitdim pressure/Θ^4 $U "Ba⋅K⁻⁴"
        @unitdim 𝟙/T/Θ $U "Hz⋅K⁻¹"
        @unitdim entropy $U "erg⋅K⁻¹"
        @unitdim specificentropy $U "erg⋅K⁻¹g⁻¹"
        @unitdim specificenergy $U "erg⋅g⁻¹"
        @unitdim thermalconductance $U "erg⋅s⁻¹K⁻¹"
        @unitdim thermalresistance $U "K⋅s⋅erg⁻¹"
        @unitdim thermalconductivity $U "erg⋅s⁻¹cm⁻¹K⁻¹"
        @unitdim thermalresistivity $U "K⋅cm⋅s⋅erg⁻¹"
    end
end

#@unitdim current EMU "Bi"
@unitdim magneticflux EMU "Mx"
@unitdim magneticfluxdensity EMU "G"
#@unitdim magneticfield EMU "Oe"
#@unitdim reluctance EMU "Bi⋅Mx⁻¹"
@unitdim magneticdipolemoment EMU "erg⋅G⁻¹"
@unitdim vectorpotential EMU "Mx⋅cm⁻¹"
#@unitdim magneticmoment EMU "Mx⋅cm"
#@unitdim polestrength EMU "pole"

#@unitdim charge Gauss "Fr"
@unitdim magneticflux Gauss "Mx"
@unitdim magneticfluxdensity Gauss "G"
#@unitdim magneticfield Gauss "Oe"
#@unitdim reluctance Gauss "Fr⋅s⁻¹Mx⁻¹"
@unitdim magneticdipolemoment Gauss "erg⋅G⁻¹"
@unitdim vectorpotential Gauss "Mx⋅cm⁻¹"
#@unitdim magneticmoment Gauss "Mx⋅cm"

@unitdim force MTS "sn"
@unitdim inv(force) MTS "sn⁻¹"
@unitdim pressure MTS "pz"
@unitdim compressibility MTS "pz⁻¹"

@unitdim mass GravitationalMetric "hyl"
@unitdim mass British "slug"
@unitdim mass IPS "slinch"
@unitdim molarmass GravitationalMetric "hyl⋅mol⁻¹"
@unitdim molarmass British "slug⋅slug-mol⁻¹"
@unitdim molarmass IPS "slinch-slinch-mol⁻¹"
@unitdim force FPS "pdl"
@unitdim pressure FPS "pdl⋅ft⁻²"
@unitdim density British "slug⋅ft⁻³"
@unitdim density IPS "slinch⋅in⁻³"
@unitdim density GravitationalMetric "hyl⋅m⁻³"

@unitdim L Rydberg "a₀"
@unitdim inv(L) Rydberg "a₀⁻¹"
@unitdim area Rydberg "a₀²"
@unitdim fuelefficiency Rydberg "a₀⁻²"
@unitdim volume Rydberg "a₀³"
@unitdim numberdensity Rydberg "a₀⁻³"
@unitdim Q Electronic "𝘦"
@unitdim Q Stoney "𝘦"
@unitdim Q Schrodinger "𝘦"
@unitdim Q CosmologicalQuantum "𝘦ₙ"
@unitdim inv(Q) Electronic "𝘦⁼¹"
@unitdim inv(Q) Stoney "𝘦⁼¹"
@unitdim inv(Q) Schrodinger "𝘦⁼¹"
@unitdim inv(Q) CosmologicalQuantum "𝘦ₙ⁼¹"
@unitdim Q^2 Electronic "𝘦²"
@unitdim Q^2 Stoney "𝘦²"
@unitdim Q^2 Schrodinger "𝘦²"
@unitdim Q^2 CosmologicalQuantum "𝘦ₙ²"
@unitdim inv(Q^2) Electronic "𝘦⁼²"
@unitdim inv(Q^2) Stoney "𝘦⁼²"
@unitdim inv(Q^2) Schrodinger "𝘦⁼²"
@unitdim inv(Q^2) CosmologicalQuantum "𝘦ₙ⁼²"

for U ∈ (:FPS,:IPS,:British,:English,:Survey)
    @eval begin
        @unitdim frequency $U "Hz"
        @unitdim frequencydrift $U "Hz⋅s⁻¹"
        @unitdim 𝟙/T/Θ $U "Hz⋅°R⁻¹"
    end
end
for U ∈ (:FPS,:British,:English,:Survey)
    @eval begin
        @unitdim luminousintensity $U "cd"
        @unitdim illuminance $U "fc"
    end
end

@doc """
$(convertext(:length,"planck(U,S)/mass(U,S)/speed(U,S)"))

Extent of one-dimensional shape or `length` (m), unit conversion factor.

```Julia
julia> L(CGS,Metric) # m⋅cm⁻¹
$(L(CGS,Metric))

julia> L(IAU,Metric) # m⋅au⁻¹
$(L(IAU,Metric))

julia> L(English,Metric) # m⋅ft⁻¹
$(L(English,Metric))

julia> L(EnglishUS,English) # ft⋅ftUS⁻¹
$(L(EnglishUS,English))

julia> L(PlanckGauss,Metric) # m⋅ℓP⁻¹
$(L(PlanckGauss,Metric))
```
""" L, ft, ftUS

@doc """
$(convertext(:time,"length(U,S)/speed(U,S)"))

Dimension along which events are ordered or `T` (s), unit conversion factor.

```Julia
julia> T(MPH,Metric) # s⋅h⁻¹
$(T(MPH,Metric))

julia> T(IAU,Metric) # s⋅D⁻¹
$(T(IAU,Metric))

julia> T(Hubble,Metric)
$(T(Hubble,Metric))
```
""" T

@doc """
    neper(U::UnitSystem) = U(𝟏,log(𝟙))

Logarithmic unit expressing the ratio of a dimensional quanty.
```Julia
julia> neper(Metric)
$(neper(Metric))

julia> exp(neper(Metric))
$(exp(neper(Metric)))
```
""" neper

@doc """
    bel(U::UnitSystem) = U(𝟏,log10(𝟙))

Logarithmic unit expressing the ratio of a dimensional quanty.
```Julia
julia> bel(Metric)
$(bel(Metric))

julia> exp10(bel(Metric))
$(exp10(bel(Metric)))
```
""" bel

@doc """
    decibel(U::UnitSystem) = U(𝟏,logdb(𝟙))

Logarithmic unit expressing the ratio of a dimensional quanty.
```Julia
julia> decibel(Metric)
$(decibel(Metric))

julia> expdb(decibel(Metric))
$(expdb(decibel(Metric)))
```
""" decibel


text:
using GeometricIntegrators.Config
using GeometricIntegrators.Integrators
using GeometricIntegrators.Solvers
using GeometricIntegrators.Tableaus
using GeometricIntegrators.TestProblems.Oscillator
using GeometricIntegrators.Utils
using Test

set_config(:nls_solver, NewtonSolver)
set_config(:jacobian_autodiff, false)

using GeometricIntegrators.TestProblems.Oscillator: Δt, nt, refx, refq, refp

ode  = oscillator_ode()
pode = oscillator_pode()
sode = oscillator_sode()
iode = oscillator_iode()
idae = oscillator_idae()
pdae = oscillator_pdae()

@testset "$(rpad("Runge-Kutta integrators",80))" begin

    @test typeof(Integrator(ode, getTableauExplicitMidpoint(), Δt)) <: IntegratorERK
    @test typeof(Integrator(ode, getTableauCrouzeix(), Δt)) <: IntegratorDIRK
    @test typeof(Integrator(ode, getTableauImplicitMidpoint(), Δt)) <: IntegratorFIRK
    @test typeof(Integrator(ode, getTableauGLRK(1), Δt)) <: IntegratorFIRK

end

@testset "$(rpad("Explicit Runge-Kutta integrators",80))" begin

    int = Integrator(ode, getTableauExplicitEuler(), Δt)
    sol = integrate(int, nt)

    @test rel_err(sol.q, refx) < 5E-2

    int = Integrator(ode, getTableauExplicitMidpoint(), Δt)
    sol = integrate(int, nt)

    @test rel_err(sol.q, refx) < 1E-3

    int = Integrator(ode, getTableauERK4(), Δt)
    sol = integrate(int, nt)

    @test rel_err(sol.q, refx) < 5E-7

end

@testset "$(rpad("Implicit Runge-Kutta integrators",80))" begin

    int = Integrator(ode, getTableauCrouzeix(), Δt)
    sol = integrate(int, nt)

    @test rel_err(sol.q, refx) < 5E-5

    int = Integrator(ode, getTableauImplicitEuler(), Δt)
    sol = integrate(int, nt)

    @test rel_err(sol.q, refx) < 5E-2

    int = Integrator(ode, getTableauImplicitMidpoint(), Δt)
    sol = integrate(int, nt)

    @test rel_err(sol.q, refx) < 5E-4

    int = Integrator(ode, getTableauGLRK(1), Δt)
    sol = integrate(int, nt)

    @test rel_err(sol.q, refx) < 5E-4

    int = Integrator(ode, getTableauGLRK(2), Δt)
    sol = integrate(int, nt)

    @test rel_err(sol.q, refx) < 1E-7

    int = Integrator(ode, getTableauGLRK(3), Δt)
    sol = integrate(int, nt)

    @test rel_err(sol.q, refx) < 1E-11

    int = Integrator(ode, getTableauGLRK(4), Δt)
    sol = integrate(int, nt)

    @test rel_err(sol.q, refx) < 1E-15

    int = Integrator(ode, getTableauGLRK(5), Δt)

    @test rel_err(sol.q, refx) < 1E-15

    int = Integrator(ode, getTableauGLRK(6), Δt)
    sol = integrate(int, nt)

    @test rel_err(sol.q, refx) < 1E-15

    int = Integrator(ode, getTableauGLRK(7), Δt)
    sol = integrate(int, nt)

    @test rel_err(sol.q, refx) < 5E-15

    int = Integrator(ode, getTableauSRK3(), Δt)
    sol = integrate(int, nt)

    @test rel_err(sol.q, refx) < 1E-7

end

@testset "$(rpad("Partitioned Runge-Kutta integrators",80))" begin

    pint = Integrator(pode, getTableauSymplecticEulerA(), Δt)
    psol = integrate(pint, nt)

    @test rel_err(psol.q, refq) < 5E-2
    @test rel_err(psol.p, refp) < 1E-3

    pint = Integrator(pode, getTableauSymplecticEulerB(), Δt)
    psol = integrate(pint, nt)

    @test rel_err(psol.q, refq) < 5E-2
    @test rel_err(psol.p, refp) < 1E-3

    pint = Integrator(pode, TableauEPRK(:prk4, 4, getTableauERK4().q), Δt)
    psol = integrate(pint, nt)

    @test rel_err(psol.q, refq) < 5E-7
    @test rel_err(psol.p, refp) < 5E-7

    pint = Integrator(pode, TableauIPRK(:pglrk, 2, getCoefficientsGLRK(1)), Δt)
    psol = integrate(pint, nt)

    @test rel_err(psol.q, refq) < 5E-4
    @test rel_err(psol.p, refp) < 5E-4

    pint = Integrator(pode, TableauIPRK(:pglrk, 4, getCoefficientsGLRK(2)), Δt)
    psol = integrate(pint, nt)

    @test rel_err(psol.q, refq) < 1E-7
    @test rel_err(psol.p, refp) < 1E-7

    pint = Integrator(pode, TableauIPRK(:pglrk, 6, getCoefficientsGLRK(3)), Δt)
    psol = integrate(pint, nt)

    @test rel_err(psol.q, refq) < 1E-11
    @test rel_err(psol.p, refp) < 1E-11

end

@testset "$(rpad("Special Runge-Kutta integrators",80))" begin

    pgint = IntegratorPGLRK(iode, getCoefficientsPGLRK(2), Δt)
    pgsol = integrate(pgint, nt)

    @test rel_err(pgsol.q, refx) < 1E-5

end

ext: jl


text:
using Rotolo

using Base.Markdown

@redirect Base.Markdown.MD

@session d1

md"test *Markdown* string"

@container bloc1 style => "font-size: small;"

md"bloc1 *Markdown* string"

@container bloc2 style => "background-color: antiquewhite;"

md"bloc2 *Markdown* string"

sleep(3)

@container bloc3 style => "background-color: lightpink;"

md"""
## title2
"""


text:
function create_logical_device(physical_device::VkPhysicalDevice, cqueue_infos::Container; features = nothing)
    device_ref = Ref{VkDevice}()
    if isnothing(features)
        features = Ref(VkPhysicalDeviceFeatures(values(DEFAULT_VK_PHYSICAL_DEVICE_FEATURES)...))
    end
    queue_infos = cqueue_infos[]
    @warn "VkDeviceCreateInfo"
    cdevice_info = @safe_contain Ref(VkDeviceCreateInfo(VK_STRUCTURE_TYPE_DEVICE_CREATE_INFO, C_NULL, 0, length(queue_infos), unsafe_pointer(queue_infos), 0, C_NULL, 0, C_NULL, unsafe_pointer(features)))
    append!(cdevice_info.deps, cqueue_infos.deps)
    @warn "vkCreateDevice"
    return
    @vkcheck vkCreateDevice(physical_device, cdevice_info[], C_NULL, device_ref)
    @warn "Device created"
    cdevice = Container(device_ref[], [cdevice_info, cqueue_infos])
    finalizer(x->vkDestroyDevice(x[], C_NULL), cdevice)
    @warn "Return"

    cdevice
end


text:
#=
# Weight normalization

The authors of <https://arxiv.org/abs/1602.07868> introduce weight normalization to boost
learning.
Let's try it here.

Preliminaries:
=#

import CairoMakie
import MLDatasets
import RestrictedBoltzmannMachines as RBMs
import WeightNormalizedRBMs as wRBMs
nothing #hide

#=
The following is a convenience function to plot grids of digits.
Given a four dimensional tensor `A` of size `(width, height, ncols, nrows)`
containing `width x height` images in a grid of `nrows x ncols`, this returns
a matrix of size `(width * ncols, height * nrows)`, that can be plotted in a heatmap
to display all images.
=#

function imggrid(A::AbstractArray{<:Any,4})
    reshape(permutedims(A, (1,3,2,4)), size(A,1)*size(A,3), size(A,2)*size(A,4))
end

#=
Load data. Keep only 0, 1 digits for speed.
=#

Float = Float32
train_x, train_y = MLDatasets.MNIST.traindata()
tests_x, tests_y = MLDatasets.MNIST.testdata()
selected_digits = (0, 1)
train_x = Array{Float}(train_x[:, :, train_y .∈ Ref(selected_digits)] .> 0.5)
tests_x = Array{Float}(tests_x[:, :, tests_y .∈ Ref(selected_digits)] .> 0.5)
train_y = train_y[train_y .∈ Ref(selected_digits)]
tests_y = tests_y[tests_y .∈ Ref(selected_digits)]
train_nsamples = length(train_y)
tests_nsamples = length(tests_y)
(train_nsamples, tests_nsamples)

#=
Init RBM and train.
=#

rbm = RBMs.BinaryRBM(zeros(Float,28,28), zeros(Float,400), zeros(Float,28,28,400))
RBMs.initialize!(rbm, train_x)
wrbm = wRBMs.WeightNormRBM(rbm) # weight normalization reparameterization
@time history = RBMs.pcd!(wrbm, train_x; epochs=500, batchsize=256, steps=5)
rbm = RBMs.RBM(wrbm)
nothing #hide

#=
Let's see what the learning curves look like.
=#

fig = CairoMakie.Figure(resolution=(800, 300))
ax = CairoMakie.Axis(fig[1,1], xlabel="train time", ylabel="log(pseudolikelihood)")
CairoMakie.lines!(ax, get(history, :lpl)...)
fig

#=
Let's look at some samples generated by this RBM.
=#

nrows, ncols = 10, 15
fantasy_x = train_x[:, :, rand(1:train_nsamples, nrows * ncols)]
@time fantasy_x .= RBMs.sample_v_from_v(rbm, fantasy_x; steps=10000)
nothing #hide

#

fig = CairoMakie.Figure(resolution=(40ncols, 40nrows))
ax = CairoMakie.Axis(fig[1,1], yreversed=true)
CairoMakie.image!(ax, imggrid(reshape(fantasy_x, 28, 28, ncols, nrows)), colorrange=(1,0))
CairoMakie.hidedecorations!(ax)
CairoMakie.hidespines!(ax)
fig


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

const Solver = Analysis
const AbstractSolver = AbstractAnalysis

#=
function Solver{S<:AbstractSolver}(::Type{S}, name="solver", properties...)
    variant = S(properties...)
    solver = Solver{S}(name, 0.0, [], [], 0, nothing, false, [], [], 0.0, Dict(), variant)
    return solver
end
=#

function Solver(::Type{S}, problems::Problem...) where S<:AbstractSolver
    solver = Solver(S, "$(S)Solver")
    push!(solver.problems, problems...)
    return solver
end

function push!(solver::Solver, problem::Problem)
    push!(solver.problems, problem)
end

function getindex(solver::Solver, problem_name::String)
    for problem in get_problems(solver)
        if problem.name == problem_name
            return problem
        end
    end
    throw(KeyError(problem_name))
end

function haskey(solver::Solver, field_name::String)
    return haskey(solver.fields, field_name)
end

get_field_problems(solver::Solver) = filter(is_field_problem, get_problems(solver))
get_boundary_problems(solver::Solver) = filter(is_boundary_problem, get_problems(solver))

"""Return one combined field assembly for a set of field problems.

Parameters
----------
solver :: Solver

Returns
-------
M, K, Kg, f, fg :: SparseMatrixCSC

Notes
-----
If several field problems exists, they are simply summed together, so
problems must have unique node ids.

"""
function get_field_assembly(solver::Solver)
    problems = get_field_problems(solver)

    problem = problems[1]
    M = problem.assembly.M
    K = problem.assembly.K
    f = problem.assembly.f
    K_csc = problem.assembly.K_csc
    f_csc = problem.assembly.f_csc
    Kg = problem.assembly.Kg
    fg = problem.assembly.fg

    for problem in problems[2:end]
        append!(M, problem.assembly.M)
        append!(K, problem.assembly.K)
        append!(Kg, problem.assembly.Kg)
        append!(f, problem.assembly.f)
        # Use in place addition with .+= ?
        K_csc += problem.assembly.K_csc
        f_csc += problem.assembly.f_csc
        append!(fg, problem.assembly.fg)
    end

    N = size(K, 1)

    M = sparse(M, N, N)
    K = sparse(K, N, N)
    if nnz(K) == 0
        @warn("Field assembly seems to be empty. Check that elements are ",
              "pushed to problem and formulation is correct.")
    end
    f = sparse(f, N, 1)
    Kg = sparse(Kg, N, N)
    fg = sparse(fg, N, 1)

    return M, problem.assemble_csc ? K_csc : K, Kg, problem.assemble_csc ? f_csc : f, fg
end

""" Loop through boundary assemblies and check for possible overconstrain situations. """
function check_for_overconstrained_dofs(solver::Solver)
    overdetermined = false
    constrained_dofs = Set{Int}()
    all_overconstrained_dofs = Set{Int}()
    boundary_problems = get_boundary_problems(solver)
    for problem in boundary_problems
        new_constraints = Set(problem.assembly.C2.I)
        new_constraints = setdiff(new_constraints, problem.assembly.removed_dofs)
        overconstrained_dofs = intersect(constrained_dofs, new_constraints)
        all_overconstrained_dofs = union(all_overconstrained_dofs, overconstrained_dofs)
        if length(overconstrained_dofs) != 0
            @warn("problem is overconstrained, finding overconstrained dofs... ")
            overdetermined = true
            for dof in overconstrained_dofs
                for problem_ in boundary_problems
                    new_constraints_ = Set(problem_.assembly.C2.I)
                    new_constraints_ = setdiff(new_constraints_, problem_.assembly.removed_dofs)
                    if dof in new_constraints_
                        @warn("overconstrained dof $dof defined in problem $(problem_.name)")
                    end
                end
                @warn("To solve overconstrained situation, remove dofs from problems so that it exists only in one.")
                @warn("To do this, use push! to add dofs to remove to problem.assembly.removed_dofs, e.g.")
                @warn("`push!(bc.assembly.removed_dofs, $dof`)")
            end
        end
        constrained_dofs = union(constrained_dofs, new_constraints)
    end
    if overdetermined
        @warn("List of all overconstrained dofs:")
        @warn(sort(collect(all_overconstrained_dofs)))
        error("problem is overconstrained, not continuing to solution.")
    end
    return true
end

""" Return one combined boundary assembly for a set of boundary problems.

Returns
-------
K, C1, C2, D, f, g :: SparseMatrixCSC

"""
function get_boundary_assembly(solver::Solver, N)

    check_for_overconstrained_dofs(solver)

    K = spzeros(N, N)
    C1 = spzeros(N, N)
    C2 = spzeros(N, N)
    D = spzeros(N, N)
    f = spzeros(N, 1)
    g = spzeros(N, 1)
    for problem in get_boundary_problems(solver)
        assembly = problem.assembly
        K_ = sparse(assembly.K, N, N)
        C1_ = sparse(assembly.C1, N, N)
        C2_ = sparse(assembly.C2, N, N)
        D_ = sparse(assembly.D, N, N)
        f_ = sparse(assembly.f, N, 1)
        g_ = sparse(assembly.g, N, 1)
        for dof in assembly.removed_dofs
            @info("$(problem.name): removing dof $dof from assembly")
            C1_[dof,:] .= 0.0
            C2_[dof,:] .= 0.0
        end
        SparseArrays.dropzeros!(C1_)
        SparseArrays.dropzeros!(C2_)

        already_constrained = get_nonzero_rows(C2)
        new_constraints = get_nonzero_rows(C2_)
        overconstrained_dofs = intersect(already_constrained, new_constraints)
        if length(overconstrained_dofs) != 0
            @warn("overconstrained dofs $overconstrained_dofs")
            @warn("already constrained = $already_constrained")
            @warn("new constraints = $new_constraints")
            overconstrained_dofs = sort(overconstrained_dofs)
            error("overconstrained dofs, not solving problem.")
        end

        K .+= K_
        C1 .+= C1_
        C2 .+= C2_
        D .+= D_
        f .+= f_
        g .+= g_
    end
    return K, C1, C2, D, f, g
end

"""
Solve linear system using LDLt factorization (SuiteSparse). This version
requires that final system is symmetric and positive definite, so boundary
conditions are first eliminated before solution.
"""
function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{1}})

    nnz(D) == 0 || return false

    A = get_nonzero_rows(K)
    B = get_nonzero_rows(C2)
    B2 = get_nonzero_columns(C2)
    B == B2 || return false
    I = setdiff(A, B)

    if length(B) == 0
        @warn("No rows in C2, forget to set Dirichlet boundary conditions to model?")
    else
        u[B] = lu(C2[B,B2]) \ Vector(g[B])
    end

    # solve interior domain using LDLt factorization
    F = ldlt(K[I,I])
    u[I] = F \ Vector(f[I] - K[I,B]*u[B])

    # solve lagrange multipliers
    la[B] = lu(C1[B2,B]) \ Vector(f[B] - K[B,I]*u[I] - K[B,B]*u[B])

    return true
end

"""
Solve linear system using LU factorization (UMFPACK). This version solves
directly the saddle point problem without elimination of boundary conditions.
It is assumed that C1 == C2 and D = 0, so problem is symmetric and zero rows
cand be removed from total system before solution. This kind of system arises
in e.g. mesh tie problem
"""
function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{2}})

    C1 == C2 || return false
    length(D) == 0 || return false

    A = [K C1'; C2  D]
    b = [f; g]
    ndofs = size(K, 2)

    nz1 = get_nonzero_rows(A)
    nz2 = get_nonzero_columns(A)
    nz1 == nz2 || return false

    x = zeros(2*ndofs)
    x[nz1] = lufact(A[nz1,nz2]) \ full(b[nz1])

    u[:] = x[1:ndofs]
    la[:] = x[ndofs+1:end]

    return true
end

"""
Solve linear system using LU factorization (UMFPACK). This version solves
directly the saddle point problem without elimination of boundary conditions.
If matrix has zero rows, diagonal term is added to that matrix is invertible.
"""
function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{3}})

    A = [K C1'; C2  D]
    b = [f; g]

    ndofs = size(K, 2)
    nonzero_rows = zeros(2*ndofs)
    for j in rowvals(A)
        nonzero_rows[j] = 1.0
    end
    A += sparse(Diagonal(1.0 .- nonzero_rows))

    x = lu(A) \ Vector(b[:])

    u[:] .= x[1:ndofs]
    la[:] .= x[ndofs+1:end]

    return true
end

""" Default linear system solver for solver. """
function solve!(solver::Solver; empty_assemblies_before_solution=true, symmetric=true)

    @info("Solving linear system.")
    t0 = Base.time()

    # assemble field & boundary problems
    # TODO: return same kind of set for both assembly types
    # M1, K1, Kg1, f1, fg1, C11, C21, D1, g1 = get_field_assembly(solver)
    # M2, K2, Kg2, f2, fg2, C12, C22, D2, g2 = get_boundary_assembly(solver)

    M, K, Kg, f, fg = get_field_assembly(solver)
    N = size(K, 2)
    Kb, C1, C2, D, fb, g = get_boundary_assembly(solver, N)
    K = K + Kg + Kb
    f = f + fg + fb

    if symmetric
        K = 1/2*(K + K')
        M = 1/2*(M + M')
    end

    if empty_assemblies_before_solution
        # free up some memory before solution by emptying field assemblies from problems
        for problem in get_field_problems(solver)
            empty!(problem.assembly)
        end
    end

    #=
    if !haskey(solver, "fint")
        solver.fields["fint"] = field(solver.time => f)
    else
        update!(solver.fields["fint"], solver.time => f)
    end

    fint = solver.fields["fint"]

    if length(fint) > 1
        # kick in generalized alpha rule for time integration
        alpha = solver.alpha
        K = (1-alpha)*K
        C1 = (1-alpha)*C1
        f = (1-alpha)*f + alpha*fint.data[end-1].second
    end
    =#

    ndofs = N
    u = zeros(ndofs)
    la = zeros(ndofs)
    is_solved = false
    local i
    for i in [1, 2, 3]
        is_solved = solve!(solver, K, C1, C2, D, f, g, u, la, Val{i})
        if is_solved
            t1 = round(Base.time()-t0; digits=2)
            norms = (norm(u), norm(la))
            @info("Solved linear system in $t1 seconds using solver $i. " *
                  "Solution norms (||u||, ||la||): $norms.")
            break
        end
    end
    if !is_solved
        error("Failed to solve linear system!")
    end
    #push!(solver.norms, norms)
    #solver.u = u
    #solver.la = la

    @info("")

    return u, la
end

"""
    assemble!(solver; with_mass_matrix=false)

Default assembler for solver.

This function loops over all problems defined in problem and launches
standard assembler for them. As a result, each problem.assembly is
populated with global stiffness matrix, force vector, and, optionally,
mass matrix.
"""
function assemble!(solver::Solver, time::Float64; with_mass_matrix=false)
    @info("Assembling problems ...")

    for problem in get_problems(solver)
        timeit("assemble $(problem.name)") do
            empty!(problem.assembly)
            assemble!(problem, time)
        end
    end

    if with_mass_matrix
        for problem in get_field_problems(solver)
            timeit("assemble $(problem.name) mass matrix") do
                assemble!(problem, time, Val{:mass_matrix})
            end
        end
    end

    #=
    ndofs = 0
    for problem in solver.problems
        Ks = size(problem.assembly.K, 2)
        Cs = size(problem.assembly.C1, 2)
        ndofs = max(ndofs, Ks, Cs)
    end
    solver.ndofs = ndofs
    =#
    @info("Assembly done!")
end

function get_unknown_fields(solver::Solver)
    fields = Dict()
    for problem in get_field_problems(solver)
        field_name = get_unknown_field_name(problem)
        field_dim = get_unknown_field_dimension(problem)
        fields[field_name] = field_dim
    end
    return fields
end

function get_unknown_field_name(solver::Solver)
    fields = get_unknown_fields(solver)
    return join(sort(collect(keys(fields))), ", ")
end

function get_unknown_field_dimension(solver::Solver)
    fields = get_unknown_fields(solver)
    return sum(values(fields))
end

""" Default initializer for solver. """
function initialize!(solver::Solver)
    if solver.initialized
        @warn("initialize!(): solver already initialized")
        return
    end
    @info("Initializing solver ...")
    problems = get_problems(solver)
    length(problems) != 0 || error("Empty solver, add problems to solver using push!")
    t0 = Base.time()
    field_problems = get_field_problems(solver)
    length(field_problems) != 0 || @warn("No field problem found from solver, add some..?")
    field_name = get_unknown_field_name(solver)
    field_dim = get_unknown_field_dimension(solver)
    @info("initialize!(): looks we are solving $field_name, $field_dim dofs/node")
    nodes = Set{Int64}()
    for problem in problems
        initialize!(problem, solver.time)
        for element in get_elements(problem)
            conn = get_connectivity(element)
            push!(nodes, conn...)
        end
    end
    nnodes = length(nodes)
    @info("Total number of nodes in problems: $nnodes")
    maxdof = maximum(nodes)*field_dim
    @info("# of max dof (=size of solution vector) is $maxdof")
    solver.u = zeros(maxdof)
    solver.la = zeros(maxdof)
    # TODO: this could be used to initialize elements too...
    # TODO: cannot initialize to zero always, construct vector from elements.
    for problem in problems
        problem.assembly.u = zeros(maxdof)
        problem.assembly.la = zeros(maxdof)
        # initialize(problem, ....)
    end
    t1 = round(Base.time()-t0; digits=2)
    @info("Initialized solver in $t1 seconds.")
    solver.initialized = true
end

function get_all_elements(solver::Solver)
    elements = [get_elements(problem) for problem in get_problems(solver)]
    return [elements...;]
end

"""
Return nodal field from all problems defined in solver.

Examples
--------
To return e.g. geometry defined in nodal points at time t=0.0, one can write:

julia> solver("geometry", 0.0)
"""
function (solver::Solver)(field_name::String, time::Float64)
    fields = []
    for problem in get_problems(solver)
        field = problem(field_name, time)
        if field == nothing
            continue
        end
        if length(field) == 0
            @warn("no field $field_name found for problem $(problem.name)")
            continue
        end
        push!(fields, field)
    end
    if length(fields) == 0
        return Dict{Integer, Vector{Float64}}()
    end
    return merge(fields...)
end

function update!(solver::Solver{S}, u, la, time) where S
    for problem in get_problems(solver)
        assembly = get_assembly(problem)
        elements = get_elements(problem)
        # update solution, first for assembly (u,la) ...
        update!(problem, assembly, u, la)
        # .. and then from assembly (u,la) to elements
        update!(problem, assembly, elements, time)
    end
end

""" Default postprocess for solver. Loop all problems and run postprocess
functions to calculate secondary fields, i.e. contact pressure, stress,
heat flux, reaction force etc. quantities.
"""
function postprocess!(solver::Solver, time)
    problems = get_problems(solver)
    nproblems = length(problems)
    @info("Postprocessing $nproblems problems.")
    for problem in problems
        for field_name in problem.postprocess_fields
            field = Val{Symbol(field_name)}
            @info("Running postprocess for problem $(problem.name), field $field_name")
            postprocess!(problem, time, field)
        end
    end
end

"""
    write_results!(solver, time)

Default xdmf update for solver. Loop all problems and write them individually
to Xdmf file. By default write the main unknown field (displacement, temperature,
...) and any fields requested separately in `problem.postprocess_fields` vector
(stress, strain, ...)
"""
function write_results!(solver, time)
    results_writers = get_results_writers(solver)
    if length(results_writers) == 0
        @info("No result writers are attached to analysis, not writing output.")
        @info("To write results to Xdmf file, attach Xdmf to analysis, i.e.")
        @info("xdmf_output = Xdmf(\"simulation_results\")")
        @info("add_results_writer!(analysis, xdmf_output)")
        return
    end
    # FIXME: result writer can be anything, not only Xdmf
    for xdmf in results_writers
        for problem in get_problems(solver)
            fields = [get_unknown_field_name(problem); problem.postprocess_fields]
            if is_boundary_problem(problem)
                fields = [fields; get_parent_field_name(problem)]
            end
            update_xdmf!(xdmf, problem, time, fields)
        end
    end
end

### Nonlinear quasistatic solver

mutable struct Nonlinear <: AbstractSolver
    time :: Float64
    iteration :: Int        # iteration counter
    min_iterations :: Int64 # minimum number of iterations
    max_iterations :: Int64 # maximum number of iterations
    convergence_tolerance :: Float64
    error_if_no_convergence :: Bool # throw error if no convergence
end

function Nonlinear()
    solver = Nonlinear(0.0, 0, 1, 10, 5.0e-5, true)
    return solver
end

""" Check convergence of problems.

Notes
-----
Default convergence criteria is obtained by checking each sub-problem convergence.
"""
function has_converged(solver::Solver{Nonlinear})
    properties = solver.properties
    converged = true
    eps = properties.convergence_tolerance
    for problem in get_field_problems(solver)
        has_converged = problem.assembly.u_norm_change < eps
        if isapprox(norm(problem.assembly.u), 0.0)
            # trivial solution
            has_converged = true
        end
        converged &= has_converged
    end
    return converged
end

""" Default solver for quasistatic nonlinear problems. """
function FEMBase.run!(solver::Solver{Nonlinear})

    time = solver.properties.time
    problems = get_problems(solver)
    properties = solver.properties

    # 1. initialize each problem so that we can start nonlinear iterations
    for problem in problems
        initialize!(problem, time)
    end

    # 2. start non-linear iterations
    for properties.iteration=1:properties.max_iterations
        @info(repeat("-", 80))
        @info("Starting nonlinear iteration #$(properties.iteration)")
        @info("Increment time t=$(round(time; digits=3))")
        @info(repeat("-", 80))

        # 2.1 update assemblies
        for problem in problems
            empty!(problem.assembly)
            assemble!(problem, time)
        end

        # 2.2 call solver for linearized system
        u, la = solve!(solver)

        # 2.3 update solution back to elements
        update!(solver, u, la, time)

        # 2.4 check convergence
        if properties.iteration >= properties.min_iterations && has_converged(solver)
            @info("Converged in $(properties.iteration) iterations.")
            # 2.4.1 run any postprocessing of problems
            postprocess!(solver, time)
            # 2.4.2 update Xdmf output
            write_results!(solver, time)
            return true
        end
    end

    # 3. did not converge
    if properties.error_if_no_convergence
        error("nonlinear iteration did not converge in $(properties.iteration) iterations!")
    end
end

### Linear quasistatic solver

""" Quasistatic solver for linear problems.

Notes
-----
Main differences in this solver, compared to nonlinear solver are:
1. system of problems is assumed to converge in one step
2. reassembly of problem is done only if it's manually requested using empty!(problem.assembly)

"""
mutable struct Linear <: AbstractSolver
    time :: Float64
end

function Linear()
    return Linear(0.0)
end

function FEMBase.run!(analysis::Analysis{Linear})
    time = analysis.properties.time
    @info("Running linear quasistatic analysis `$(analysis.name)` at time $time.")
    problems = get_problems(analysis)
    nproblems = length(problems)
    @info("Assembling $nproblems problems.")
    @timeit "assemble problems" for problem in problems
        isempty(problem.assembly) || continue
        initialize!(problem, time)
        assemble!(problem, time)
    end
    @timeit "solve linear system" u, la = solve!(analysis)
    @timeit "update problems" update!(analysis, u, la, time)
    postprocess!(analysis, time)
    write_results!(analysis, time)
    @info("Quasistatic linear analysis ready.")
end

# Convenience functions

function LinearSolver(name::String="Linear solver")
    return Solver(Linear, name)
end

function LinearSolver(problems::Problem...)
    solver = LinearSolver()
    add_problems!(solver, collect(problems))
    return solver
end

function NonlinearSolver(name::String="Nonlinear solver")
    return Solver(Nonlinear, name)
end

function NonlinearSolver(problems::Problem...)
    solver = NonlinearSolver()
    add_problems!(solver, collect(problems))
    return solver
end

# will be deprecated

function (solver::Solver)(time::Float64=0.0)
    @warn("analysis(time) is deprecated. Instead, use run!(analysis)")
    solver.properties.time = time
    run!(solver)
end

function solve!(solver::Solver, time::Float64)
    @warn("solve!(analysis, time) is deprecated. Instead, use run!(analysis)")
    solver.properties.time = time
    run!(solver)
end


text:
# ====== New stuff

# ===================
# Tuples
# ===================

const AbstractTuple = Union{Tuple, NamedTuple}

const NamedNTuple{N, T} = NamedTuple{S, NTuple{N, T}} where S
const AbstractNTuple{N, T} = Union{NTuple{N, T}, NamedNTuple{N, T}}

const TupleOfTypes = Tuple{Vararg{Type}}
const NamedTupleOfTypes = NamedTuple{S, <:Tuple{Vararg{Type}}} where {S}
const AbstractTupleOfTypes = Union{TupleOfTypes, NamedTupleOfTypes}

# TODO maybe better name
function tuple_of_types(TT::Type{<:Tuple})
    return Tuple(TT.types)
end

function tuple_of_types(TT::Type{<:NamedTuple})
    return NamedTuple{Tuple(TT.names)}(TT.types)
end

# TODO maybe better name
function to_tuple(types::TupleOfTypes)
    return Tuple{types...}
end

function typeTuple(types::NamedTupleOfTypes)
    return NamedTuple{keys(types), Tuple{values(types)...}}
end

const Adjlist{T} = Vector{Vector{T}}

# TODO better name, do we still need these?
function E_VAL_for_key(E_VALS::Type{<:Tuple}, key)
    return E_VALS.types[key]
end

function E_VAL_for_key(E_VALS::Type{<:NamedTuple}, key)
    return NamedTuple{Tuple(E_VALS.names)}(E_VALS.types)[key]
end

# ==========

function Adjlist{T}(n::Integer) where {T}
    result = Adjlist{T}(undef, n)
    @inbounds for i in 1:n
        result[i] = T[]
    end
    return result
end

function deepcopy_adjlist(adjlist::Adjlist{T}) where {T}
    n = length(adjlist)
    result = Vector{Vector{T}}(undef, n)
    @inbounds for i in OneTo(n)
        list = adjlist[i]
        n_list = length(list)
        result_list = Vector{T}(undef, n_list)
        for j in OneTo(n_list)
            result_list[j] = deepcopy(list[j])
        end
        result[i] = result_list
    end
    return result
end

#=
const EdgeValContainer{T} = Union{Adjlist{T},
                                  Tuple{Vararg{Adjlist}},
                                  NamedTuple{S, <: Tuple{Vararg{Adjlist}}} where S
                                 }
=#

# TODO remove these two?
create_edgeval_list(nv, E_VAL::TupleOfTypes) = Tuple(Adjlist{T}(nv) for T in E_VAL.types)
create_edgeval_list(nv, E_VAL::NamedTupleOfTypes) = NamedTuple{Tuple(E_VAL.names)}(Adjlist{T}(nv) for T in E_VAL.types)

# TODO check if still correct
edgevaluetype_from_edgevaluecontainertype(::Type{<:Adjlist{T}}) where {T} = T
edgevaluetype_from_edgevaluecontainertype(T::Type{<:Tuple})  =
    Tuple{ (edgevaluetype_from_edgevaluecontainertype(TT) for TT in T.types)...  }
edgevaluetype_from_edgevaluecontainertype(T::Type{<:NamedTuple})  =
    NamedTuple{ Tuple(T.names), Tuple{ (edgevaluetype_from_edgevaluecontainertype(TT) for TT in T.types)... } }

# TODO check if still correct
edgevals_container_type(::Val{Nothing}) where {E_VAL <: Type} = Adjlist{Nothing}
edgevals_container_type(::Val{E_VAL}) where {E_VAL <: Type} = Adjlist{E_VAL}

# TODO check if still correct
@generated function edgevals_container_type(::Val{E_VAL}) where {E_VAL <:Tuple}
    R = Tuple{( Adjlist{T} for T in E_VAL.types )...}
    return :($R)
end

# TODO check if still correct
@generated function edgevals_container_type(::Val{E_VAL}) where {E_VAL <:NamedTuple}
    R = NamedTuple{ Tuple(E_VAL.names), Tuple{( Adjlist{T} for T in E_VAL.types )...}}
    return :($R)
end

#= TODO remove?
function set_value_for_index!(adjlist::Adjlist,
                              s::Integer,
                              index::Integer,
                              value)
    @inbounds adjlist[s][index] = value
    return nothing
end
=#

@generated function set_values_for_index!(edgevals::AbstractTuple, s::Integer, index::Integer, values::T) where {T <: AbstractTuple}
    len = length(T.types)
    exprs = Expr[]
    for i in OneTo(len)
        e = :(@inbounds edgevals[$i][s][index] = values[$i])
        push!(exprs, e)
    end
    return Expr(:block, exprs...)
end

#=
function set_value_for_index!(tup_adjlist::AbstractTuple, s::Integer, index::Integer, key, value)
    @inbounds tup_adjlist[key][s][index] = value
    return nothing
end
=#

#= TODO remove?
function insert_value_for_index!(adjlist::Adjlist,
                                 s::Integer,
                                 index::Integer,
                                 value)
    @inbounds insert!(adjlist[s], index, value)
    return nothing
end
=#

@generated function insert_values_for_index!(tup_adjlist::AbstractTuple,
                                 s::Integer,
                                 index::Integer,
                                 values::T) where {T <: AbstractTuple}
    len = length(T.types)
    exprs = Expr[]
    for i in OneTo(len)
        e = :(@inbounds insert!(tup_adjlist[$i][s], index, values[$i]))
        push!(exprs, e)
    end
    return Expr(:block, exprs...)
end

#= TODO remove?
function delete_value_for_index!(adjlist::Adjlist,
                                 s::Integer,
                                 index::Integer)
    @inbounds deleteat!(adjlist[s], index)
    return nothing
end
=#

function delete_values_for_index!(tup_adjlist::AbstractTuple,
                                 s::Integer,
                                 index::Integer)
    @inbounds for adjlist in tup_adjlist
        deleteat!(adjlist[s], index)
    end
end

function value_for_index(adjlist::Adjlist, s::Integer, index::Integer)
    @inbounds return adjlist[s][index]
end

# TODO E_VAL could probably be extracted without being explicitly passed to the function
function values_for_index(tup_adjlist::AbstractTuple, E_VAL::Type, s::Integer, index::Integer)
    @inbounds return E_VAL( adjlist[s][index] for adjlist in tup_adjlist )
end

function values_for_index(tup_adjlist::AbstractTuple, E_VAL, s::Integer, index::Integer, key)
    adjlist = tup_adjlist[key]
    @inbounds return adjlist[s][index]
end

struct NoKey end
const nokey = NoKey()


text:
# Note: poor performance

# this may not be working properly

cd(@__DIR__)
using AdvancedPS, Random, StatsBase, Distributions

N = 30
# Parameters for each cluster, we assume that each cluster is Gaussian distributed in the example.
μs = [-3.5, 0.0]
# Construct the data points.
data = mapreduce(c -> rand(MvNormal([μs[c], μs[c]], 1.), N), hcat, 1:2)

priors = (
    idx=(Categorical([.5,.5]),N),
    μs = (Normal(0,1),2),
)

bounds = ((1,2),(0.0,Inf))

function loglike(idx, μs, data)
    LL = 0.0
    for (i,id) in enumerate(idx)
        LL += logpdf(MvNormal([μs[id], μs[id]], 1.0), data[:,i])
    end
    return LL
end

loglike(θ) = loglike(θ..., data)

model = DEModel(priors=priors, model=loglike)

de = DE(bounds=bounds, burnin=1000, priors=priors, progress=true, Np=15)
n_iter = 2000
@elapsed chains = sample(model, de, n_iter)
println(chains)

max_stars_repo_path: DEMCMC_Examples/Discrete Example.jl

max_issues_repo_path: DEMCMC_Examples/Discrete Example.jl

max_forks_repo_path: DEMCMC_Examples/Discrete Example.jl


text:
function plotdist(points, ego, idx)
    points = points'
    p1 = plot()
    scatter!(p1, points[:,1], points[:,2])
    scatter!(p1, [ego[1]], [ego[2]])
    plot!(p1, [ego[1], points[idx,:][1]], [ego[2], points[idx,:][2]])
    return p1
end

max_stars_repo_name: Affos/Bamberg

max_issues_repo_name: Affos/Bamberg

max_forks_repo_name: Affos/Bamberg


text:
using FileIO
using Printf
include("config.jl")
include("calc_loader.jl")
include("data_loader.jl")
include("utils/utils.jl")

function save_result(X, Xr, runtime, method, reqrank, cnt, datasetname)
    tensor_size = size(X)
    LSerror = norm(X-Xr)
    if method in non_negative_methods
        KLerror = get_dkl(X, Xr)
    else
        KLerror = missing
    end
    result = Dict(
                  "tensor_size" => tensor_size,
                  "runtime" => runtime,
                  "LSerror" => LSerror,
                  "KLerror" => KLerror,
                  "method" => method,
                  "reqrank" => reqrank
                 )
    reqrank_str =  get_tensorinfo(reqrank)
    savepath = "../result/real/$datasetname/reqrank$reqrank_str"
    mkpath(savepath)
    save("$savepath/$method$cnt.jld2", result)
    return LSerror, KLerror
end

function run(datasetname, reqranks)
    X = load_tensor(datasetname)
    tensor_size = size(X)
    println("input tensor $datasetname shape $tensor_size")
    println("__________________________________________")
    @printf("%5s %10s %20s %8s %10s %10s %10s \n", "cnt", "datasetname", "reqrank", "method", "runtime", "LSerror", "KLerror")
    println("__________________________________________")
    for reqrank in reqranks
        for method in methods
            for cnt = 1:rep_max
                runtime = @elapsed begin
                    Xr = calc(X, method, reqrank)
                end
                LSerror, KLerror = save_result(X, Xr, runtime, method, reqrank, cnt, datasetname)
                @printf("%5s %10s %20s %8s %10f %10f %10s \n", cnt, datasetname, reqrank, method, runtime, LSerror, KLerror)
            end
        end
        println("__________________________________________")
    end
end

datasetname = "AttFace"
run(datasetname, plans[datasetname])
#datasetname = "4DLFD"
#run(datasetname, plans[datasetname])

max_stars_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_issues_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_forks_repo_name: gkazunii/Legendre-tucker-rank-reduction


text:
input = joinpath(@__DIR__, "input")
input = read(input, String)

input = split(input, "\n\n")

one = input[1]
two = input[2]

CARDS = r"(\d+)"

one = map(m -> parse(Int, m.captures[1]), collect(eachmatch(CARDS, one)))
two = map(m -> parse(Int, m.captures[1]), collect(eachmatch(CARDS, two)))

# remove player number from matches
one = one[2:end]
two = two[2:end]

function play(one, two)
    while length(one) > 0 && length(two) > 0
        x = popfirst!(one)
        y = popfirst!(two)
        if x > y
            push!(one, x, y)
        else
            push!(two, y, x)
        end
    end
    if length(one) > length(two)
        return (one, nothing)
    else
        return (nothing, two)
    end
end

function playr(one, two)
    seen = Set{}()
    while length(one) > 0 && length(two) > 0
        if (one, two) in seen
            return (one, nothing)
        end
        push!(seen, (copy(one), copy(two)))
        x = popfirst!(one)
        y = popfirst!(two)
        if x <= length(one) && y <= length(two)
            result = playr(copy(one[1:x]), copy(two[1:y]))
            if isnothing(result[1])
                push!(two, y, x)
            else
                push!(one, x, y)
            end
        elseif x > y
            push!(one, x, y)
        else
            push!(two, y, x)
        end
    end
    if length(one) > length(two)
        return (one, nothing)
    else
        return (nothing, two)
    end
end

function score(deck)
    reverse!(deck)
    sum = 0
    for i in 1:length(deck)
        sum += i * deck[i]
    end
    return sum
end

result = play(copy(one), copy(two))
winner = isnothing(result[1]) ? result[2] : result[1]
p1 = score(winner)

result = playr(copy(one), copy(two))
winner = isnothing(result[1]) ? result[2] : result[1]
p2 = score(winner)

println("-----------------------------------------------------------------------")
println("crab combat -- part one :: $p1")
println("crab combat -- part two :: $p2")
println("-----------------------------------------------------------------------")

@assert(p1 == 32448)
@assert(p2 == 32949)


text:
export Plane

mutable struct Plane{T<:AbstractFloat} <: Robot
  mass::T
  g::T
  ρ::T
  Area::T
  Cd0::T
  Kd::T
  α_0::T
  v_hat::T    # z_dot \in [-0.1,0.1]*v_hat
  v_min::T    # v \in [3,10] m/s
  v_max::T
  γ_max::T    # gamma \in [-π/6,π/6]
  ϕ_max::T    # phi \in [-π/4,π/4]
  acc_lim::T  # u_a \in [-8,8] m/s
  ϕ_d_lim::T  # ϕ_d \in [-120,120] deg/s
  α_d_lim::T  # α_d \in [-60,60] deg/s
  btCollisionObject
end

function Plane{T}() where T
  # X: (x,y,z,psi,v,gamma,phi, alpha)
  # pos- xyz, heading- psi, speed- v, flight path angle- gamma, roll- phi, angle of attack- alpha,
  # U: (u_acc, u_ϕ_d, u_α_d)

  mass = 1.
  g = 9.81
  ρ = 1.225
  Area = 0.7
  Cd0 = 0.015
  Kd = 0.025
  α_0 = 5*π/180
  γ_max = π/6
  ϕ_max = π/4
  v_min = 3.    # min speed
  v_max = 10.   # max speed

  C_con = π*ρ*Area/mass
  v_hat = sqrt(g/(C_con*α_0))
  acc_nom = (ρ*Area)*(v_hat^2)*(Cd0+4*(π^2)*Kd*(α_0^2))

  acc_lim = 10*acc_nom      # linear acceleration lim
  ϕ_d_lim = (π/3)*2      # ϕ_dot lim
  α_d_lim = 2*γ_max # α_dot lim

  btCollisionObject = BulletCollision.sphere(SVector{3}(zeros(3)), 1.)

  return Plane{T}(mass,g,ρ,Area,Cd0,Kd,α_0,v_hat,v_min,v_max,γ_max,ϕ_max,acc_lim,ϕ_d_lim,α_d_lim,btCollisionObject)
end
Plane(::Type{T} = Float64; kwargs...) where {T} = Plane{T}(; kwargs...)

function configuration(x, rb::Plane)
  # convert state to (x,y,z,yaw,pitch,roll)
  pitch = rb.α_0 - x[8] - x[6]
  return [x[1:4]; pitch; x[7]]
end

function body2inertial(x,rb)
  # See 2.2 in Small Unmanned Aircraft: Theory and Practice by Beard & McLain

  X = configuration(x,rb) # X: (x,y,z,yaw,pitch,roll)

  R_z = [cos(yaw) -sin(yaw) 0;
         sin(yaw) cos(yaw) 0;
         0 0 1]

  pitch = -pitch  # since not using NED

  R_y = [cos(pitch) 0 sin(pitch);
         0 1 0;
        -sin(pitch) 0 cos(pitch)]

  R_x = [1 0 0;
         0 cos(roll) -sin(roll);
         0 sin(roll) cos(roll)]

  R = R_z*R_y*R_x   # R: body->inertial
end

function inertial2body(x,rb)
  return body2inertial(x,rb)' # R: inertial->body
end

function plane_geometry()
  α_0 = 5*π/180
  λ = 0.1

  # Vertices
  body_V =
      [10. -0.5 -1.
       10.  0.5 -1.
       10. -0.5  0.
       10.  0.5  0.
        5. -0.5  1.
        5.  0.5  1.
      -10. -0.5  1.
      -10.  0.5  1.
      -10. -0.5  0.
      -10.  0.5  0.
      -7.5 -0.5 -1.
      -7.5  0.5 -1.]*λ
  body = FlexibleConvexHull([Vector(body_V[k,:]) for k in 1:size(body_V,1)])

  wing_V = [ 3.   0.   3*sin(α_0)+.1
             3.   0.   3*sin(α_0)-.1
             0. -10.  .1
             0. -10. -.1
             0.  10.  .1
             0.  10. -.1
            -2. -10.  -2*sin(α_0)+.1
            -2. -10.  -2*sin(α_0)-.1
            -2.  10.  -2*sin(α_0)+.1
            -2.  10.  -2*sin(α_0)-.1]*λ
  wing = FlexibleConvexHull([Vector(wing_V[k,:]) for k in 1:size(wing_V,1)])

  stab_V = [ -8.  0. .5
             -8.  0. .1
             -9. -4. .5
             -9. -4. .1
             -9.  4. .5
             -9.  4. .1
             -10. -4. .5
             -10. -4. .1
             -10.  4. .5
             -10.  4. .1]*λ
  stab = FlexibleConvexHull([Vector(stab_V[k,:]) for k in 1:size(stab_V,1)])

  tail_V =  [-7. -.2 1.
             -6. .2 1.
              -10. -.2 1.
              -10. .2 1.
              -10. -.2 4.
              -10. .2 4.
              -9. -.2 4.
              -9. .2 4.]*λ
  tail = FlexibleConvexHull([Vector(tail_V[k,:]) for k in 1:size(tail_V,1)])

  box_V = [-1 -1 -0.1
          1 -1 -0.1
          1 1 -0.1
          -1 1 -0.1
          -1 -1 0.4
          1 -1 0.4
          1  1 0.4
          -1 1 0.4]
  box = FlexibleConvexHull([Vector(box_V[k,:]) for k in 1:size(box_V,1)])

  return body, wing, stab, tail, box
end


text:
#!/opt/local/bin/julia

# include("bin/mean.jl")

using DelimitedFiles
using Statistics
using StatsBase
using Distributions
using JSON
using Plots

data = readdlm("sample.ssv", ' ');
rates = Array{Float64}( data[:,5] )

µ = Statistics.mean( rates )
σ = Statistics.std( rates )
γ = StatsBase.skewness( rates, µ )
κ = StatsBase.kurtosis( rates, µ )

z = StatsBase.zscore( rates, µ, σ )
y = map( x -> round(x,digits=1), z )

println( "sample µ $µ  σ $σ  γ $γ  κ $κ" )

dto = Dict(
       "mean" => µ,
       "std" => σ,
       "data" => z,
      )

open( "out.json", "w" ) do io
    write( io, JSON.json(dto) )
end

d = Distributions.fit( Normal, rates )
println( d )

# Plots.histogram( z )

# h = StatsBase.fit( Histogram, z )
# Plots.plot( h )

max_stars_repo_name: karlredgate/oaf

max_issues_repo_name: karlredgate/oaf

max_forks_repo_name: karlredgate/oaf


text:
##############################################################################
#
# Specialization of AbstractMPCCModel when mp, G, H are given by NLPModels
#
# TODO:
# - weird cons function
# - jac, jacG, jacH
# - check hess operators
#
##############################################################################

"""
Definit le type MPCC :
min f(x)
l <= x <= u
lb <= c(x) <= ub
lccG <= G(x) _|_ H(x) >= lccH
"""
mutable struct MPCCNLPs <: AbstractMPCCModel

 mp :: AbstractNLPModel
 G  :: AbstractNLPModel
 H  :: AbstractNLPModel

 meta :: MPCCModelMeta

 counters :: MPCCCounters

 function MPCCNLPs(mp  :: AbstractNLPModel;
                   G   :: AbstractNLPModel = ADNLPModel(x->0, mp.meta.x0),
                   H   :: AbstractNLPModel = ADNLPModel(x->0, mp.meta.x0),
                   ncc :: Int64 = -1,
                   x0  :: Vector = mp.meta.x0)

  ncc = ncc == -1 ? G.meta.ncon : ncc

  if G.meta.ncon != H.meta.ncon || (ncc != -1 && G.meta.ncon != ncc)
   throw(error("Incompatible complementarity"))
  end

  n    = length(x0)
  ncon = mp.meta.ncon

  meta = MPCCModelMeta(n, x0 = x0, ncc = ncc, lccG = G.meta.lcon, lccH = H.meta.lcon,
                              lvar = mp.meta.lvar, uvar = mp.meta.uvar,
                              ncon = ncon,
                              lcon = mp.meta.lcon, ucon = mp.meta.ucon)

  return new(mp, G, H, meta, MPCCCounters())
 end
end

############################################################################

# Getteur

############################################################################
"""
Evaluate f(x), the objective function at x.
"""
function obj(mod :: MPCCNLPs, x :: AbstractVector)
 increment!(mod, :neval_obj)
 return obj(mod.mp, x)
end

"""
Evaluate ∇f(x), the gradient of the objective function at x in place.
"""
function grad!(mod :: MPCCNLPs, x :: Vector, gx :: AbstractVector)
 increment!(mod, :neval_grad)
 return grad!(mod.mp, x, gx)#gx
end

"""
Evaluate ``c(x)``, the constraints at `x`.
"""
function cons_nl!(mod :: MPCCNLPs, x :: Vector, c :: AbstractVector)
 increment!(mod, :neval_cons)
 return cons!(mod.mp, x, c)
end

"""
Evaluate ``G(x)``, the constraints at `x`.
"""
function consG!(mod :: MPCCNLPs, x :: Vector, c :: AbstractVector)
 increment!(mod, :neval_consG)
 if mod.meta.ncc > 0
  c .= cons(mod.G, x)
 else
  c .= Float64[]
 end

 return c
end

"""
Evaluate ``H(x)``, the constraints at `x`.
"""
function consH!(mod :: MPCCNLPs, x :: Vector, c :: AbstractVector)
 increment!(mod, :neval_consH)
 if mod.meta.ncc > 0
  c .= cons(mod.H, x)
 else
  c .= Float64[]
 end

 return c
end

"""
    Jx = jac(nlp, x)
Evaluate ``[∇c(x), -∇G(x), -∇H(x)]``, the constraint's Jacobian at `x` as a sparse matrix.
Size of the jacobian matrix: (ncon + 2 * ncc) x n
"""
function jac(mod :: MPCCNLPs, x :: AbstractVector)
  n, ncon, ncc = mod.meta.nvar, mod.meta.ncon, mod.meta.ncc

 if ncon+ncc == 0

  A = sparse(zeros(0,n))

 else

  if ncon > 0
   J = jac_nl(mod, x)
  else
   J = sparse(zeros(0,n))
  end

  if ncc > 0
   JG, JH = jacG(mod,x), jacH(mod,x)
   A = [J;-JG;-JH]
  else
   A = J
  end

 end

 return A
end

"""
Evaluate ``∇c(x)``, the constraint's Jacobian at `x` as a sparse matrix.
"""
function jac_nl(mod :: MPCCNLPs, x :: AbstractVector)
 increment!(mod, :neval_jac)
 return jac(mod.mp, x)
end

"""
Evaluate ``∇G(x)``, the constraint's Jacobian at `x` as a sparse matrix.
"""
function jacG(mod :: MPCCNLPs, x :: AbstractVector)
 increment!(mod, :neval_jacG)
 if mod.meta.ncc > 0
  rslt = jac(mod.G, x)
 else
  rslt = Float64[]
 end

 return rslt
end

"""
Evaluate ``∇H(x)``, the constraint's Jacobian at `x` as a sparse matrix.
"""
function jacH(mod :: MPCCNLPs, x :: AbstractVector)
 increment!(mod, :neval_jacH)
 if mod.meta.ncc > 0
  rslt = jac(mod.H, x)
 else
  rslt = Float64[]
 end

 return rslt
end

"""
Jacobienne des contraintes actives à precmpcc près

(Tangi: copie de colonne si contrainte d'égalité?)
"""
function jacl(mod :: MPCCNLPs, x :: Vector) #bad name
 A, Il,Iu,Ig,Ih,IG,IH = jac_actif(mod, x, 0.0)
 return A
end

function jac_actif(mod :: MPCCNLPs, x :: Vector, prec :: Float64)

  n = mod.meta.nvar
  ncc = mod.meta.ncc

  Il = findall(z->z<=prec,abs.(x-mod.meta.lvar))
  Iu = findall(z->z<=prec,abs.(x-mod.meta.uvar))
  jl = zeros(n);jl[Il] .= 1.0;Jl=diagm(0 => jl);
  ju = zeros(n);ju[Iu] .= 1.0;Ju=diagm(0 => ju);
  IG=Int64[];IH=Int64[];Ig=Int64[];Ih=Int64[];

 if mod.meta.ncon+ncc ==0

  A=[]

 else

  if mod.meta.ncon > 0
   c = cons_nl(mod,x)
   J = jac_nl(mod, x)
  else
   c = Float64[]
   J = sparse(zeros(0,2))
  end

  Ig=findall(z->z<=prec,abs.(c-mod.meta.lcon))
  Ih=findall(z->z<=prec,abs.(c-mod.meta.ucon))

  Jg, Jh = zeros(mod.meta.ncon,n), zeros(mod.meta.ncon,n)

  Jg[Ig,1:n] = J[Ig,1:n]
  Jh[Ih,1:n] = J[Ih,1:n]

  if ncc>0

   IG=findall(z->z<=prec,abs.(consG(mod,x)-mod.meta.lccG))
   IH=findall(z->z<=prec,abs.(consH(mod,x)-mod.meta.lccH))

   JG, JH = zeros(ncc, n), zeros(ncc, n)
   JG[IG,1:n] = jacG(mod,x)[IG,1:n]
   JH[IH,1:n] = jacH(mod,x)[IH,1:n]

   A=[Jl;Ju;-Jg;Jh;-JG;-JH]'
  else
   A=[Jl;Ju;-Jg;Jh]'
  end

 end

 return A, Il,Iu,Ig,Ih,IG,IH
end

"""
    Jv = jnlprod(nlp, x, v, Jtv)
Evaluate ``∇c(x)v``, the transposed-Jacobian-vector product at `x`.
"""
function jnlprod!(mod :: MPCCNLPs, x :: Vector, v :: Vector, Jv :: Vector)
 return jprod!(mod.mp, x, v, Jv)
end

"""
  JGv = jGprod!(nlp, x, v, Jv)
Evaluate ``∇G(x)v``, the Jacobian-vector product at `x` in place.
"""
function jGprod!(mod :: MPCCNLPs, x :: AbstractVector, v :: AbstractVector, Jv :: AbstractVector)
    increment!(mod, :neval_jGprod)
    if mod.meta.ncc > 0
     Jv .= jprod(mod.G, x, v)
    else
     Jv .= Float64[]
    end
    return Jv
end

"""
  JHv = jHprod!(nlp, x, v, Jv)
Evaluate ``∇H(x)v``, the Jacobian-vector product at `x` in place.
"""
function jHprod!(mod :: MPCCNLPs, x :: AbstractVector, v :: AbstractVector, Jv :: AbstractVector)
    increment!(mod, :neval_jHprod)
    if mod.meta.ncc > 0
     Jv .= jprod(mod.H, x, v)
    else
     Jv .= Float64[]
    end
    return Jv
end

"""
    Jtv = jtprod(nlp, x, v, Jtv)
Evaluate ``∇c(x)^Tv``, the transposed-Jacobian-vector product at `x`.
"""
function jnltprod!(mod :: MPCCNLPs, x :: Vector, v :: Vector, Jv :: Vector)
 increment!(mod, :neval_jtprod)
 return jtprod!(mod.mp, x, v, Jv)
end

"""
    Jtv = jtprodG(nlp, x, v, Jtv)
Evaluate ``∇G(x)^Tv``, the transposed-Jacobian-vector product at `x`
"""
function jGtprod!(mod :: MPCCNLPs, x :: Vector, v :: Vector, Jv :: Vector)
 increment!(mod, :neval_jGtprod)
 if mod.meta.ncc > 0
  Jv .= jtprod(mod.G, x, v)
 else
  Jv = Float64[]
 end

 return Jv
end

"""
    Jtv = jtprodH(nlp, x, v, Jtv)
Evaluate ``∇H(x)^Tv``, the transposed-Jacobian-vector product at `x`
"""
function jHtprod!(mod :: MPCCNLPs, x :: Vector, v :: Vector, Jv :: Vector)
 increment!(mod, :neval_jHtprod)
 if mod.meta.ncc > 0
  Jv .= jtprod(mod.H, x, v)
 else
  Jv = Float64[]
 end

 return Jv
end

function hess(mod :: MPCCNLPs, x :: Vector ; obj_weight = 1.0)
 increment!(mod, :neval_hess)
 return hess(mod.mp, x, obj_weight = obj_weight)
end

function hess(mod :: MPCCNLPs, x :: AbstractVector, y :: AbstractVector; obj_weight = 1.0)
 increment!(mod, :neval_hess)
 ncon, ncc = mod.meta.ncon, mod.meta.ncc
 yG, yH = y[ncon+1:ncon+ncc], y[ncon+ncc+1:ncon+2*ncc]
 return hess(mod.mp, x, y[1:ncon], obj_weight = obj_weight) - hessG(mod, x, yG) - hessH(mod, x, yH)
end

function hessG(mod :: MPCCNLPs, x :: AbstractVector, y :: AbstractVector; obj_weight = 1.0)
 increment!(mod, :neval_hessG)
 if mod.meta.ncc > 0
  rslt = hess(mod.G, x, y, obj_weight = obj_weight)
 else
  rslt = zeros(0,0)
 end
 return rslt
end

function hessH(mod :: MPCCNLPs, x :: AbstractVector, y :: AbstractVector; obj_weight = 1.0)
 increment!(mod, :neval_hessH)
 if mod.meta.ncc > 0
  rslt = hess(mod.H, x, y, obj_weight = obj_weight)
 else
  rslt = zeros(0,0)
 end
 return rslt
end

function hess_structure!(mod :: MPCCNLPs, rows :: AbstractVector{<: Integer}, cols :: AbstractVector{<: Integer})
  n = mod.meta.nvar
  I = ((i,j) for i = 1:n, j = 1:n if i ≥ j)
  rows .= getindex.(I, 1)
  cols .= getindex.(I, 2)
  return rows, cols
end

function hess_coord!(nlp :: MPCCNLPs, x :: AbstractVector, vals :: AbstractVector; obj_weight :: Real = one(eltype(x)))
  increment!(nlp, :neval_hess)
  Hx = hess(nlp, x, obj_weight=obj_weight)
  k = 1
  for j = 1 : nlp.meta.nvar
    for i = j : nlp.meta.nvar
      vals[k] = Hx[i, j]
      k += 1
    end
  end
  return vals
end

function hess_coord!(nlp :: MPCCNLPs, x :: AbstractVector, y :: AbstractVector, vals :: AbstractVector; obj_weight :: Real = one(eltype(x)))
  increment!(nlp, :neval_hess)
  Hx = hess(nlp, x, y, obj_weight=obj_weight)
  k = 1
  for j = 1 : nlp.meta.nvar
    for i = j : nlp.meta.nvar
      vals[k] = Hx[i, j]
      k += 1
    end
  end
  return vals
end

function hprod!(nlp :: MPCCNLPs, x :: AbstractVector, v :: AbstractVector, Hv :: AbstractVector; obj_weight :: Real = one(eltype(x)))
  increment!(nlp, :neval_hprod)
  Hx = hess(nlp, x, obj_weight = obj_weight)
  Hv .= (Hx + Hx' - diagm(0 => diag(Hx))) * v
  return Hv
end

function hprod!(nlp :: MPCCNLPs, x :: AbstractVector, y :: AbstractVector, v :: AbstractVector, Hv :: AbstractVector; obj_weight :: Real = one(eltype(x)))
  increment!(nlp, :neval_hprod)
  Hx = hess(nlp, x, y, obj_weight = obj_weight)
  Hv .= (Hx + Hx' - diagm(0 => diag(Hx))) * v
  return Hv
end


text:
using CutPruners

@testset "DeMatos cut pruning" begin
    algo = DeMatosPruningAlgo(3)
    pruner = CutPruner{2, Int}(algo, :Max)
    addposition!(pruner, [2 0])
    addcuts!(pruner, [1 0], [0], [true])
    addcuts!(pruner, [2 0], [0], [true])
    addcuts!(pruner, [3 0], [0], [true])
    addusage!(pruner, [1 0])
    # FIXME: need to update stats before reaching the ncuts limit ...
    addposition!(pruner, [1 0])
    @test pruner.trust == [0, 0, 2.0]
    addcuts!(pruner, [4 0], [0], [true])
    @test pruner.A[:,1] == [4, 2, 3]
    @test pruner.trust == [2.0, 0, 0]
    addposition!(pruner, [1, 0])
    @test pruner.trust == [3.0, 0, 0]
    addcuts!(pruner, [5 0], [0], [true])
    @test pruner.A[:,1] == [4, 5, 3]
    @test pruner.trust == [0, 3.0, 0]

    @test ncuts(pruner) == 3

    addposition!(pruner, [1 0])
    @test pruner.trust == [0, 4.0, 0]
end

@testset "Lazy Minus" begin
    algo = DeMatosPruningAlgo(2)
    pruner = CutPruner{2, Int}(algo, :Max, true)
    addcuts!(pruner, [1 0], [0], [true])
    addcuts!(pruner, [2 0], [0], [true])
    addposition!(pruner, [1 0])
    @test pruner.trust == [1., 0.]
    addcuts!(pruner, [-1 0], [0], [true])
    @test pruner.trust == [0., 1.]
    @test pruner.A == [1 0; -1 0]
end

@testset "Priority to new cuts" begin
    algo = DeMatosPruningAlgo(1)
    pruner = CutPruner{2, Int}(algo, :Min)
    addcuts!(pruner, [1 0], [0], [true])
    addposition!(pruner, [1 0])
    addcuts!(pruner, [2 0], [0], [true])

    @test ncuts(pruner) == 1
    # The cut [2 0] is taken even if [1 0] has one territory
    @test pruner.A == [2 0]
end


text:
macro import_canvas_huge()
    :(
        using ComplexVisual:
            CV_Canvas, CV_2DCanvas, cv_create_context, cv_save_image,
            cv_destroy, CV_Std2DCanvas, CV_Math2DCanvas,
            cv_pixel2math, cv_math2pixel
    )
end

import Base:show

"""
`CV_Canvas`: something where you can draw and paint.
"""
abstract type CV_Canvas  end

"""
```
CV_2DCanvas <: CV_Canvas
    # required fields:
    # surface       (CairoSurface)
    # pixel_width   (convert(Int32, surface.width))
    # pixel_height  (convert(Int32, surface.height))
    # bounding_box  (CV_Rectangle{Int32})  bottom=0, left=0 (i.e. zero-based)
```

a thin wrapper around a `CairoSurface` (where `pixel_width` and `pixel_height`
are `Int32`, instead of `Floats`) for twodimensional painting and drawing
operations.

There are canvases with "only" pixel coorindates (like `CV_Std2DCanvas`)
and there are canvases which also have a mathematical coordinate
system (like `CV_Math2DCanvas`).
"""
abstract type CV_2DCanvas <: CV_Canvas    # {{{
    # required fields:
    # surface       (CairoSurface)
    # pixel_width   (convert(Int32, surface.width))
    # pixel_height  (convert(Int32, surface.height))
    # bounding_box  (CV_Rectangle{Int32})  bottom=0, left=0 (i.e. zero-based)
end

show(io::IO, canvas::CV_2DCanvas) = cv_show_impl(io, canvas)
show(io::IO, m::MIME{Symbol("text/plain")}, canvas::CV_2DCanvas) =
    cv_show_impl(io, m, canvas)

function cv_destroy(can::CV_2DCanvas)
    destroy(can.surface)
    return nothing
end

function cv_anchor(can::CV_2DCanvas, anchor_name::Symbol)
    return cv_anchor(can.bounding_box, anchor_name)
end

"""
```
cv_create_context(canvas; prepare=true)
    canvas    CV_2DCanvas
    prepare   Bool
```

Create Context for 2DCanvas. If `prepare` is `true` then
(depending on the concrete CV_2DCanvas subtype) special preperations
or initializations for the context are done.
"""
function cv_create_context(canvas::CV_2DCanvas; prepare::Bool=true)
    return CV_2DCanvasContext(canvas)
end

# }}}

"""
```
cv_create_cairo_image_surface(width, height)
    width     Integer
    height    Integer
```
Use ARGB32 image surface with pixel matrix stored in julia.
"""
function cv_create_cairo_image_surface(width::Integer, height::Integer)
    return CairoImageSurface(
        zeros(UInt32, width, height),
        Cairo.FORMAT_ARGB32; flipxy=false)
end

"""
```
cv_save_image(canvas, filename)
    canvas      CV_2DCanvas
    filename    AbstractString
```

save content of canvas as a png image.
"""
function cv_save_image(canvas::CV_2DCanvas, filename::AbstractString)
    write_to_png(canvas.surface, filename)
    return nothing
end

"""
```
CV_Std2DCanvas <: CV_2DCanvas
    surface          Cairo.CairoSurfaceImage{UInt32}
    pixel_width      Int32
    pixel_height     Int32
    bounding_box     CV_Rectangle{Int32}
```

A twodimensional canvas with a `CairoSurfaceImage` (ARGB32 format).
"""
struct CV_Std2DCanvas <: CV_2DCanvas  # {{{
    surface      :: Cairo.CairoSurfaceImage{UInt32}
    pixel_width  :: Int32
    pixel_height :: Int32
    bounding_box :: CV_Rectangle{Int32}

    function CV_Std2DCanvas(pixel_width::Integer, pixel_height::Integer)
        width = convert(Int32, pixel_width)
        height = convert(Int32, pixel_height)
        surface = cv_create_cairo_image_surface(width, height)
        self = new(
            surface, width, height,
            CV_Rectangle(height, Int32(0), Int32(0), width))
        return self
    end
end # }}}

"""
```
CV_Math2DCanvas <: CV_2DCanvas
    corner_ul       Complex{Float64}  # (math.) coordinates of upper left corner
    corner_lr       Complex{Float64}  # (math.) coordinates of lower right corner
    resolution      Float64           # Pixels per (math.) unit
    surface         Cairo.CairoSurfaceImage{UInt32}
    pixel_width     Int32
    pixel_height    Int32
    bounding_box    CV_Rectangle{Int32}
```

A twodimensional canvas with a `CairoSurfaceImage` (ARGB32 format) where
the user coordinates represent a mathmatical coordinate system.
With the two methods `cv_math2pixel` and `cv_pixel2math` one can
transform math coordinates to pixel coordinates and vice versa.
See also: `CV_MathCoorStyle`.
"""
struct CV_Math2DCanvas <: CV_2DCanvas # {{{
    corner_ul    :: Complex{Float64}; # (math.) coordinates of upper left corner
    corner_lr    :: Complex{Float64}; # (math.) coordinates of lower right corner
    resolution   :: Float64;          # Pixels per (math.) unit
    surface      :: Cairo.CairoSurfaceImage{UInt32}
    pixel_width  :: Int32
    pixel_height :: Int32
    bounding_box :: CV_Rectangle{Int32}

    function CV_Math2DCanvas(
            corner_ul::Complex{Float64}, corner_lr::Complex{Float64},
            resolution=100.0::Float64)    # {{{
        if (!isfinite(resolution))
            cv_error("resolution must be finite; resolution = ", resolution)
        end
        if (!(resolution > 0))
            cv_error("resolution must be positive; resolution = ", resolution)
        end
        if (!isfinite(corner_ul) || !isfinite(corner_lr))
            cv_error(
                "corner_ul and corner_lr must be finite",
                "; corner_ul = ", corner_ul,
                "; corner_lr = ", corner_lr)
        end
        if (!(real(corner_ul) < real(corner_lr)))
            cv_error(
                "requirement: real(corner_ul) < real(corner_lr)",
                "; corner_ul = ", corner_ul,
                "; corner_lr = ", corner_lr)
        end
        if (!(imag(corner_ul) > imag(corner_lr)))
            cv_error(
                "requirement: imag(corner_ul) > imag(corner_lr)",
                "; corner_ul = ", corner_ul,
                "; corner_lr = ", corner_lr)
        end
        math_width = real(corner_lr) - real(corner_ul)
        math_height = imag(corner_ul) - imag(corner_lr)
        pixel_width = floor(Int32, math_width*resolution)
        pixel_height = floor(Int32, math_height*resolution)
        if (pixel_width == 0 || pixel_height == 0)
            cv_error(
                "resolution too small; image width or height would be 0",
                "; math_width = ", math_width,
                "; math_height = ", math_height,
                "; pixel_width = ", pixel_width,
                "; pixel_height = ", pixel_height)
        end
        surface = CairoImageSurface(
            zeros(UInt32, pixel_width, pixel_height),
            Cairo.FORMAT_ARGB32; flipxy=false)
        self = new(
            corner_ul, corner_lr, resolution, surface,
            pixel_width, pixel_height,
            CV_Rectangle(pixel_height, Int32(0), Int32(0), pixel_width))
        return self
    end  # }}}
end

function show(io::IO, canvas::CV_Math2DCanvas)
    t = typeof(canvas)
    print(io, string(t.name.name), "(⌜")
    show(io, canvas.corner_ul)
    print(io, ", ")
    show(io, canvas.corner_lr)
    print(io, "⌟, res=")
    show(io, canvas.resolution)
    print(io, ", ")
    show(io, canvas.bounding_box)
    print(io, ')')
    return nothing
end

show(io::IO, m::MIME{Symbol("text/plain")}, canvas::CV_Math2DCanvas) =
    show(io, canvas)

"""
```
cv_math2pixel(canvas, mx, my) :: Tuple{Int32, Int32}
    canvas      CV_Math2DCanvas
    mx          Float64
    my          Float64
```

convert math coordinates `(mx, my)` to pixel coordinates.
"""
function cv_math2pixel(canvas::CV_Math2DCanvas, mx::Float64, my::Float64)
    res = canvas.resolution
    return (
        round(Int32, (mx - real(canvas.corner_ul))*res),
        round(Int32, -(my - imag(canvas.corner_ul))*res))
end

"""
```
cv_pixel2math(canvas, px, py::Integer) :: Tuple{Float64, Float64}
    canvas    CV_Math2DCanvas
    px        Integer
    pyi       Integer
```

convert pixel coordinates `(px, py)` to math coordinates.
"""
function cv_pixel2math(canvas::CV_Math2DCanvas, px::Integer, py::Integer)
    res = canvas.resolution
    return (
        real(canvas.corner_ul) + px/res,
        imag(canvas.corner_ul) - py/res)
end

"""
```
cv_create_context(canvas; prepare=true)
    canvas       CV_Math2DCanvas
    prepare      Bool
```
create context for Math2DCanvas. If `prepare` is `true` the
the context's user coordinate system is the mathematical coordinate
system of the Math2DCanvas. See also `CV_MathCoorStyle`.
"""
function cv_create_context(canvas::CV_Math2DCanvas; prepare::Bool=true)
    con = CV_2DCanvasContext(canvas)
    if prepare
        cv_prepare(con, CV_MathCoorStyle(canvas))
    end
    return con
end

# }}}

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
export CIF

mutable struct CIF
    volume::Float64
    a::Float64
    b::Float64
    c::Float64
end

mutable struct CIF
    comments::Vector{String}

end

function readcif(filename::String)
    lines = open(readlines, filename)
    for line in lines
        if line[1] == '#' # comment

        end
        items = split(line, " ", keepempty=false)
        if startswith(line, "data_")
            data = line[6:end]
        elseif items[1] == "loop_"

        elseif items[1] == "_cell_length_a"
            x.a = parse(Float64, items[2])
        elseif items[1] == "_cell_length_b"
            x.b = parse(Float64, items[2])
        elseif items[1] == "_cell_length_c"
            x.c = parse(Float64, items[2])
        elseif items[1] == "_cell_angle_alpha"
            x.α = parse(Float64, items[2])
        elseif items[1] == "_cell_angle_beta"
            x.β = parse(Float64, items[2])
        elseif items[1] == "_cell_angle_gamma"
            x.γ = parse(Float64, items[2])
        elseif items[1] == "_cell_volume"
            x.volume = parse(Float64, items[2])
        elseif "# Dielectric constant, electronic:"
        else

        end
    end
end

function parseaaa(lines::Vector{String})
    for i = 1:7
        @assert startswith(lines[i], " _atom_site")
    end
    for i = 8:length(lines)
        line = strip(lines[i])
        isempty(line) && continue
        items = split(line, " ", keepempty=false)
        atom = items[1]
        x = parse(Float64, items[4])
        y = parse(Float64, items[5])
        z = parse(Float64, items[6])

    end
end