text:
module DimensionalAnalysis

export PhysicalQuantity
export Mass, Length, Time, Temperature, Intensity, Current, Amount
export Dimensionless
export Area, Volume, Density, Speed, Acceleration, Force, Energy, Power, Stress, Pressure
export Diffusivity, DynamicViscosity, KinematicViscosity, SurfaceEnergy

export mass, length_, time_, temperature, intensity, current, amount
export dimensionless
export area, volume, density, speed, acceleration, force, energy, power, stress, pressure
export diffusivity, kinematicViscosity, dynamicViscosity, surfaceEnergy

export ReynoldsNumber, WeberNumber, ModifiedWeberNumber

export MeasuredQuantity
export Meter

type PhysicalQuantity{X<:Real,Mass,Length,Time,Temperature,Intensity,Current,Amount}
    value::X
end

value{X,Mass,Length,Time,Temperature,Intensity,Current,Amount}(x::PhysicalQuantity{X,Mass,Length,Time,Temperature,Intensity,Current,Amount} ) = x.value

+{X,Y,Mass,Length,Time,Temperature,Intensity,Current,Amount}(x::PhysicalQuantity{X,Mass,Length,Time,Temperature,Intensity,Current,Amount},
                                                             y::PhysicalQuantity{Y,Mass,Length,Time,Temperature,Intensity,Current,Amount}) =
    PhysicalQuantity{typeof(x.value+y.value),Mass,Length,Time,Temperature,Intensity,Current,Amount}(x.value+y.value)

-{X,Y,Mass,Length,Time,Temperature,Intensity,Current,Amount}(x::PhysicalQuantity{X,Mass,Length,Time,Temperature,Intensity,Current,Amount},
                                                             y::PhysicalQuantity{Y,Mass,Length,Time,Temperature,Intensity,Current}) =
    PhysicalQuantity{typeof(x.value-y.value),Mass,Length,Time,Temperature,Intensity,Current,Amount}(x.value-y.value)

*{X,MassX,LengthX,TimeX,TemperatureX,IntensityX,CurrentX,AmountX,Y,MassY,LengthY,TimeY,TemperatureY,IntensityY,CurrentY,AmountY}(x::PhysicalQuantity{X,MassX,LengthX,TimeX,TemperatureX,IntensityX,CurrentX,AmountX},
                                                                                y::PhysicalQuantity{Y,MassY,LengthY,TimeY,TemperatureY,IntensityY,CurrentY,AmountY}) =
    PhysicalQuantity{typeof(x.value*y.value),MassX+MassY,LengthX+LengthY,TimeX+TimeY,TemperatureX+TemperatureY,IntensityX+IntensityY,CurrentX+CurrentY,AmountX+AmountY}(x.value*y.value)

/{X,MassX,LengthX,TimeX,TemperatureX,IntensityX,CurrentX,AmountX,Y,MassY,LengthY,TimeY,TemperatureY,IntensityY,CurrentY,AmountY}(x::PhysicalQuantity{X,MassX,LengthX,TimeX,TemperatureX,IntensityX,CurrentX,AmountX},
                                                                                y::PhysicalQuantity{Y,MassY,LengthY,TimeY,TemperatureY,IntensityY,CurrentY,AmountY}) =
    PhysicalQuantity{typeof(x.value/y.value),MassX-MassY,LengthX-LengthY,TimeX-TimeY,TemperatureX-TemperatureY,IntensityX-IntensityY,CurrentX-CurrentY,AmountX-AmountY}(x.value/y.value)

typealias Dimensionless{X} PhysicalQuantity{X, 0, 0, 0, 0, 0, 0, 0}

dimensionless{X}( x::X ) = Dimensionless{X}( x )
dimensionless() = dimensionless( 1 )

<{X<:Number,Y}(x::X, y::Dimensionless{Y})  = x < value( y )
<{X,Y<:Number}(x::Dimensionless{X}, y::Y)   = value( x ) < y

<={X<:Number,Y}(x::X, y::Dimensionless{Y}) = x <= value( y )
<={X,Y<:Number}(x::Dimensionless{X}, y::Y)  = value( x ) <= y

=={X<:Number,Y}(x::X, y::Dimensionless{Y}) = x == value( y )
=={X,Y<:Number}(x::Dimensionless{X}, y::Y)  = value( x ) == y

!={X<:Number,Y}(x::X, y::Dimensionless{Y}) = x != value( y )
!={X,Y<:Number}(x::Dimensionless{X}, y::Y) = value( x ) != y

>={X<:Number,Y}(x::X, y::Dimensionless{Y}) = x >= value( y )
>={X,Y<:Number}(x::Dimensionless{X}, y::Y)  = value( x ) >= y

>{X<:Number,Y}(x::X, y::Dimensionless{Y})  = x > value( y )
>{X,Y<:Number}(x::Dimensionless{X}, y::Y)   = value( x ) > y

*{X,Mass,Length,Time,Temperature,Intensity,Current,Amount,Y<:Number}( x::PhysicalQuantity{X,Mass,Length,Time,Temperature,Intensity,Current,Amount}, y::Y ) = x*dimensionless( y )
*{X<:Number,Y,Mass,Length,Time,Temperature,Intensity,Current,Amount}( x::X, y::PhysicalQuantity{Y,Mass,Length,Time,Temperature,Intensity,Current,Amount} ) = dimensionless( x )*y

/{X,Mass,Length,Time,Temperature,Intensity,Current,Amount,Y}( x::PhysicalQuantity{X,Mass,Length,Time,Temperature,Intensity,Current,Amount}, y::Y ) = x/dimensionless( y )
/{X, Y,Mass,Length,Time,Temperature,Intensity,Current,Amount}( x::X, y::PhysicalQuantity{Y,Mass,Length,Time,Temperature,Intensity,Current,Amount} ) = dimensionless( x )/y

>{X,Y,Mass,Length,Time,Temperature,Intensity,Current,Amount}( x::PhysicalQuantity{X,Mass,Length,Time,Temperature,Intensity,Current,Amount}, y::PhysicalQuantity{Y,Mass,Length,Time,Temperature,Intensity,Current,Amount} ) = value( x ) >  value( y )
>={X,Y,Mass,Length,Time,Temperature,Intensity,Current,Amount}( x::PhysicalQuantity{X,Mass,Length,Time,Temperature,Intensity,Current,Amount}, y::PhysicalQuantity{Y,Mass,Length,Time,Temperature,Intensity,Current,Amount}) = value( x ) >= value( y )
=={X,Y,Mass,Length,Time,Temperature,Intensity,Current,Amount}( x::PhysicalQuantity{X,Mass,Length,Time,Temperature,Intensity,Current,Amount}, y::PhysicalQuantity{Y,Mass,Length,Time,Temperature,Intensity,Current,Amount}) = value( x ) == value( y )
!={X,Y,Mass,Length,Time,Temperature,Intensity,Current,Amount}( x::PhysicalQuantity{X,Mass,Length,Time,Temperature,Intensity,Current,Amount}, y::PhysicalQuantity{Y,Mass,Length,Time,Temperature,Intensity,Current,Amount}) = value( x ) != value( y )
<={X,Y,Mass,Length,Time,Temperature,Intensity,Current,Amount}( x::PhysicalQuantity{X,Mass,Length,Time,Temperature,Intensity,Current,Amount}, y::PhysicalQuantity{Y,Mass,Length,Time,Temperature,Intensity,Current,Amount}) = value( x ) <= value( y )
<{X,Y,Mass,Length,Time,Temperature,Intensity,Current,Amount}( x::PhysicalQuantity{X,Mass,Length,Time,Temperature,Intensity,Current,Amount}, y::PhysicalQuantity{Y,Mass,Length,Time,Temperature,Intensity,Current,Amount} ) = value( x ) <  value( y )

typealias Mass{X} PhysicalQuantity{X, 1, 0, 0, 0, 0, 0, 0}
mass{X}( x::X ) = Mass{X}( x )
mass() = mass( 1 )

typealias Length{X} PhysicalQuantity{X, 0, 1, 0, 0, 0, 0, 0}
length_{X}( x::X ) = Length{X}( x )
length_() = length_( 1 )

typealias Time{X} PhysicalQuantity{X, 0, 0, 1, 0, 0, 0, 0}
time_{X}( x::X ) = Time{X}( x )
time_() = time_( 1 )

typealias Temperature{X} PhysicalQuantity{X, 0, 0, 0, 1, 0, 0, 0}
temperature{X}( x::X ) = Temperature{X}( x )
temperature() = temperature( 1 )

typealias Intensity{X} PhysicalQuantity{X, 0, 0, 0, 0, 1, 0, 0}
intensity{X}( x::X ) = Intensity{X}( x  )
intensity() = intensity( 1 )

typealias Current{X} PhysicalQuantity{X, 0, 0, 0, 0, 0, 1, 0}
current{X}( x::X ) = Current{X}( x )
current() = current( 1 )

typealias Amount{X} PhysicalQuantity{X, 0, 0, 0, 0, 0, 0, 1}
amount{X}( x::X ) = Amount{X}( x )
amount() = amount( 1 )

typealias Area{X} PhysicalQuantity{X, 0, 2, 0, 0, 0, 0, 0}
area{X}( x::X ) = Area{X}( x )
area() = area( 1 )

typealias Volume{X} PhysicalQuantity{X, 0, 3, 0, 0, 0, 0, 0}
volume{X}( x::X ) = Volume{X}( x )
volume() = volume( 1 )

typealias Density{X} PhysicalQuantity{X, 1,-3, 0, 0, 0, 0, 0}
density{X}( x::X ) = Density{X}( x )
density() = density( 1 )

typealias Speed{X} PhysicalQuantity{X, 0, 1,-1, 0, 0, 0, 0}
speed{X}( x::X ) = Speed{X}( x )
speed() = speed( 1 )

typealias Acceleration{X} PhysicalQuantity{X, 0, 1,-2, 0, 0, 0, 0}
acceleration{X}( x::X ) = Acceleration{X}( x )
acceleration() = acceleration( 1 )

typealias Force{X} PhysicalQuantity{X, 1, 1,-2, 0, 0, 0, 0}
force{X}( x::X ) = Force{X}( x )
force() = force( 1 )

typealias Momentum{X} PhysicalQuantity{X, 1, 1,-1, 0, 0, 0, 0}
momentum{X}( x::X ) = Momentum{X}( x )
momentum() = momentum( 1 )

typealias Energy{X} PhysicalQuantity{X, 1, 2,-2, 0, 0, 0, 0}
energy{X}( x::X ) = Energy{X}( x )
energy() = energy( 1 )

typealias Power{X} PhysicalQuantity{X, 1, 2,-3, 0, 0, 0, 0}
power{X}( x::X ) = Power{X}( x )
power() = power( 1 )

typealias Stress{X} PhysicalQuantity{X, 1,-1,-2, 0, 0, 0, 0}
stress{X}( x::X ) = Stress{X}( x )
stress() = stress( 1 )

typealias Pressure{X} PhysicalQuantity{X, 1,-1,-2, 0, 0, 0, 0}
pressure{X}( x::X ) = Pressure{X}( x )
pressure() = pressure( 1 )

typealias Diffusivity{X} PhysicalQuantity{X, 0, 2,-1, 0, 0, 0, 0}
diffusivity{X}( x::X ) = Diffusivity{X}( x )
diffusivity() = diffusivity( 1 )

typealias KinematicViscosity{X} PhysicalQuantity{X, 0, 2, -1, 0, 0, 0, 0}
kinematicViscosity{X}( x::X ) = KinematicViscosity{X}( x )
kinematicViscosity() = kinematicViscosity( 1 )

typealias DynamicViscosity{X} PhysicalQuantity{X, 1,-1,-1, 0, 0, 0, 0}
dynamicViscosity{X}( x::X ) = DynamicViscosity{X}( x )
dynamicViscosity() = dynamicViscosity( 1 )

typealias SurfaceEnergy{X} PhysicalQuantity{X, 1, 0,-2, 0, 0, 0, 0}
surfaceEnergy{X}( x::X ) = SurfaceEnergy{X}( x )
surfaceEnergy() = surfaceEnergy( 1 )

ReynoldsNumber{X1,X2,X3,X4}( ρ::Density{X1}, v::Speed{X2}, d::Length{X3}, μ::DynamicViscosity{X4} ) = ρ*v*d/μ
ReynoldsNumber{X1,X2,X3}( v::Speed{X1}, d::Length{X2}, ν::KinematicViscosity{X3 } ) = v*d/ν

WeberNumber{X1,X2,X3,X4}( ρ::Density{X1}, v::Speed{X2}, d::Length{X3}, σ::SurfaceEnergy{X4} ) = ρ*v^2*d/σ

ModifiedWeberNumber{X1,X2,X3,X4}( ρ::Density{X1}, v::Speed{X2}, d::Length{X3}, σ::SurfaceEnergy{X4} ) = WeberNumber( ρ, v, d, σ )/12

type MeasuredQuantity{T,Mass,Length,Time,Temperature,Intensity,Current,Amount,Offset,Scale}
    value::PhysicalQuantity{T,Mass,Length,Time,Temperature,Intensity,Current,Amount}
end

value{T,Mass,Length,Time,Temperature,Intensity,Current,Amount,Offset,Scale}( x::MeasuredQuantity{T,Mass,Length,Time,Temperature,Intensity,Current,Amount,Offset,Scale}) =
    value( x.value )

typealias Meter{T} MeasuredQuantity{T,0,1,0,0,0,0,0,0,1}

show{T}(   x::Meter{T}) = println( value( x ), "m" )
offset{T}( x::Meter{T}) = 0
scale{T}(  x::Meter{T}) = 1

# M,L,T,Θ,I,C,N

# SchmidtNumber{T1,T2,T3}( ρ::MechanicalDensity{T1}, μ::MechanicalDynamicViscosity{T2}, d::MechanicalMolecularDiffusivity{T}) = μ*d/ρ
# SchmidtNumber{T1,T2}( ν::MechanicalKinematicViscosity{T1}, d::MechanicalMolecularDiffusivity{T2} ) = ν/d

# # Weber number	We	\mathrm{We} = \frac{\rho v^2 l}{\sigma}	multiphase flow (strongly curved surfaces; ratio of inertia to surface tension)
# WeberNumber{T1,T2,T3,T4}( ρ::Mechanical

# Abbe number	V	V = \frac{ n_d - 1 }{ n_F - n_C }	optics (dispersion in optical materials)
# Activity coefficient	\gamma	 \gamma= \frac {{a}}{{x}} 	chemistry (Proportion of "active" molecules or atoms)
# Albedo	\alpha	\alpha= (1-D) \bar \alpha(\theta_i) + D \bar{ \bar \alpha}	climatology, astronomy (reflectivity of surfaces or bodies)
# Archimedes number	Ar	 \mathrm{Ar} = \frac{g L^3 \rho_\ell (\rho - \rho_\ell)}{\mu^2}	fluid mechanics (motion of fluids due to density differences)
# Arrhenius number	\alpha	\alpha = \frac{E_a}{RT} 	chemistry (ratio of activation energy to thermal energy)[5]
# Atomic weight	M		chemistry (mass of atom over one atomic mass unit, u, where carbon-12 is exactly 12 u)
# Atwood number	A	\mathrm{A} = \frac{\rho_1 - \rho_2} {\rho_1 + \rho_2} 	fluid mechanics (onset of instabilities in fluid mixtures due to density differences)
# Bagnold number	Ba	\mathrm{Ba} = \frac{\rho d^2 \lambda^{1/2} \gamma}{\mu}	fluid mechanics, geology (ratio of grain collision stresses to viscous fluid stresses in flow of a granular material such as grain and sand)[6]
# Bejan number

# (fluid mechanics)	Be	\mathrm{Be} = \frac{\Delta P L^2} {\mu \alpha}	fluid mechanics (dimensionless pressure drop along a channel)[7]
# Bejan number
# (thermodynamics)	Be	\mathrm{Be} = \frac{\dot S'_{\mathrm{gen},\, \Delta T}}{\dot S'_{\mathrm{gen},\, \Delta T}+ \dot S'_{\mathrm{gen},\, \Delta p}}	thermodynamics (ratio of heat transfer irreversibility to total irreversibility due to heat transfer and fluid friction)[8]
# Bingham number	Bm	\mathrm{Bm} = \frac{ \tau_y L }{ \mu V }	fluid mechanics, rheology (ratio of yield stress to viscous stress)[5]
# Biot number	Bi	\mathrm{Bi} = \frac{h L_C}{k_b}	heat transfer (surface vs. volume conductivity of solids)
# Blake number	Bl or B	\mathrm{B} = \frac{u \rho}{\mu (1 - \epsilon) D}	geology, fluid mechanics, porous media (inertial over viscous forces in fluid flow through porous media)
# Bodenstein number	Bo or Bd	\mathrm{Bo} = vL/\mathcal{D} = \mathrm{Re}\, \mathrm{Sc} 	chemistry (residence-time distribution; similar to the axial mass transfer Peclet number)[9]
# Bond number	Bo	\mathrm{Bo} = \frac{\rho a L^2}{\gamma}	geology, fluid mechanics, porous media (buoyant versus capilary forces, similar to the Eötvös number) [10]
# Brinkman number	Br	 \mathrm{Br} = \frac {\mu U^2}{\kappa (T_w - T_0)}	heat transfer, fluid mechanics (conduction from a wall to a viscous fluid)
# Brownell–Katz number	NBK	\mathrm{N}_\mathrm{BK} = \frac{u \mu}{k_\mathrm{rw}\sigma} 	fluid mechanics (combination of capillary number and Bond number) [11]
# Capillary number	Ca	\mathrm{Ca} = \frac{\mu V}{\gamma} 	porous media, fluid mechanics (viscous forces versus surface tension)
# Chandrasekhar number	Q	 \mathrm{Q} = \frac{{B_0}^2 d^2}{\mu_0 \rho \nu \lambda} 	magnetohydrodynamics (ratio of the Lorentz force to the viscosity in magnetic convection)
# Colburn J factors	JM, JH, JD		turbulence; heat, mass, and momentum transfer (dimensionless transfer coefficients)
# Coefficient of kinetic friction	\mu_k		mechanics (friction of solid bodies in translational motion)
# Coefficient of static friction	\mu_s		mechanics (friction of solid bodies at rest)
# Coefficient of determination	R^2		statistics (proportion of variance explained by a statistical model)
# Coefficient of variation	\frac{\sigma}{\mu}	\frac{\sigma}{\mu}	statistics (ratio of standard deviation to expectation)
# Correlation	ρ or r	\frac{{\mathbb E}[(X-\mu_X)(Y-\mu_Y)]}{\sigma_X \sigma_Y} or  \frac{\sum_{k=1}^n (x_k-\bar x)(y_k-\bar y)}{\sqrt{\sum_{k=1}^n (x_k-\bar x)^2 \sum_{k=1}^n (y_k-\bar y)^2}}  where \bar x = \sum_{k=1}^n x_k/n and similarly for \bar y	statistics (measure of linear dependence)
# Courant–Friedrich–Levy number	C or 𝜈	C = \frac {u\,\Delta t} {\Delta x}	mathematics (numerical solutions of hyperbolic PDEs)[12]
# Damkohler number	Da	 \mathrm{Da} = k \tau	chemistry (reaction time scales vs. residence time)
# Damping ratio	\zeta	 \zeta = \frac{c}{2 \sqrt{km}}	mechanics (the level of damping in a system)
# Darcy friction factor	Cf or fD		fluid mechanics (fraction of pressure losses due to friction in a pipe; four times the Fanning friction factor)
# Darcy number	Da	 \mathrm{Da} = \frac{K}{d^2}	porous media (ratio of permeability to cross-sectional area)
# Dean number	D	\mathrm{D} = \frac{\rho V d}{\mu} \left( \frac{d}{2 R} \right)^{1/2}	turbulent flow (vortices in curved ducts)
# Deborah number	De	 \mathrm{De} = \frac{t_\mathrm{c}}{t_\mathrm{p}}	rheology (viscoelastic fluids)
# Decibel	dB		acoustics, electronics, control theory (ratio of two intensities or powers of a wave)
# Drag coefficient	cd	c_\mathrm{d} = \dfrac{2 F_\mathrm{d}}{\rho v^2 A}\, ,	aeronautics, fluid dynamics (resistance to fluid motion)
# Dukhin number	Du	 \mathrm{Du} = \frac{\kappa^{\sigma}}{{\Kappa_m} a}	colloid science (ratio of electric surface conductivity to the electric bulk conductivity in heterogeneous systems)
# Eckert number	Ec	 \mathrm{Ec} = \frac{V^2}{c_p\Delta T}  	convective heat transfer (characterizes dissipation of energy; ratio of kinetic energy to enthalpy)
# Ekman number	Ek	\mathrm{Ek} = \frac{\nu}{2D^2\Omega\sin\varphi} 	geophysics (viscous versus Coriolis forces)
# Elasticity
# (economics)	E	E_{x,y} = \frac{\partial \ln(x)}{\partial \ln(y)} = \frac{\partial x}{\partial y}\frac{y}{x}	economics (response of demand or supply to price changes)
# Eötvös number	Eo	\mathrm{Eo}=\frac{\Delta\rho \,g \,L^2}{\sigma}	fluid mechanics (shape of bubbles or drops)
# Ericksen number	Er	\mathrm{Er}=\frac{\mu v L}{K}	fluid dynamics (liquid crystal flow behavior; viscous over elastic forces)
# Euler number	Eu	 \mathrm{Eu}=\frac{\Delta{}p}{\rho V^2} 	hydrodynamics (stream pressure versus inertia forces)
# Euler's number	e	e =  \displaystyle\sum\limits_{n = 0}^{ \infty} \dfrac{1}{n!} \approx 2.71828 	mathematics (base of the natural logarithm)
# Excess temperature coefficient	\Theta_r	\Theta_r = \frac{c_p (T-T_e)}{U_e^2/2}	heat transfer, fluid dynamics (change in internal energy versus kinetic energy)[13]
# Fanning friction factor	f		fluid mechanics (fraction of pressure losses due to friction in a pipe; 1/4th the Darcy friction factor)[14]
# Feigenbaum constants	\alpha, \delta	\alpha \approx 2.50290,
# \ \delta \approx 4.66920 	chaos theory (period doubling)[15]
# Fine structure constant	\alpha	\alpha = \frac{e^2}{2\varepsilon_0 hc}	quantum electrodynamics (QED) (coupling constant characterizing the strength of the electromagnetic interaction)
# f-number	f	 f = \frac {{\ell}}{{D}}	optics, photography (ratio of focal length to diameter of aperture)
# Föppl–von Kármán number	\gamma	\gamma = \frac{Y r^2}{\kappa}	virology, solid mechanics (thin-shell buckling)
# Fourier number	Fo	\mathrm{Fo} = \frac{\alpha t}{L^2}	heat transfer, mass transfer (ratio of diffusive rate versus storage rate)
# Fresnel number	F	\mathit{F} = \frac{a^{2}}{L \lambda}	optics (slit diffraction)[16]
# Froude number	Fr	\mathrm{Fr} = \frac{v}{\sqrt{g\ell}}	fluid mechanics (wave and surface behaviour; ratio of a body's inertia to gravitational forces)
# Gain	–		electronics (signal output to signal input)
# Gain ratio	–		bicycling (system of representing gearing; length traveled over length pedaled)[17]
# Galilei number	Ga	\mathrm{Ga} = \frac{g\, L^3}{\nu^2}	fluid mechanics (gravitational over viscous forces)
# Golden ratio	\varphi	\varphi = \frac{1+\sqrt{5}}{2} \approx 1.61803	mathematics, aesthetics (long side length of self-similar rectangle)
# Görtler number	G	\mathrm{G} = \frac{U_e \theta}{\nu} \left( \frac{\theta}{R} \right)^{1/2}	fluid dynamics (boundary layer flow along a concave wall)
# Graetz number	Gz	\mathrm{Gz} = {D_H \over L} \mathrm{Re}\, \mathrm{Pr}	heat transfer, fluid mechanics (laminar flow through a conduit; also used in mass transfer)
# Grashof number	Gr	 \mathrm{Gr}_L = \frac{g \beta (T_s - T_\infty ) L^3}{\nu ^2}	heat transfer, natural convection (ratio of the buoyancy to viscous force)
# Gravitational coupling constant	\alpha_G	\alpha_G=\frac{Gm_e^2}{\hbar c}	gravitation (attraction between two massy elementary particles; analogous to the Fine structure constant)
# Hatta number	Ha	\mathrm{Ha} = \frac{N_{\mathrm{A}0}}{N_{\mathrm{A}0}^{\mathrm{phys}}} 	chemical engineering (adsorption enhancement due to chemical reaction)
# Hagen number	Hg	 \mathrm{Hg} = -\frac{1}{\rho}\frac{\mathrm{d} p}{\mathrm{d} x}\frac{L^3}{\nu^2} 	heat transfer (ratio of the buoyancy to viscous force in forced convection)
# Hydraulic gradient	i	i = \frac{\mathrm{d}h}{\mathrm{d}l} = \frac{h_2 - h_1}{\mathrm{length}}	fluid mechanics, groundwater flow (pressure head over distance)
# Iribarren number	Ir	\mathrm{Ir} = \frac{\tan \alpha}{\sqrt{H/L_0}}	wave mechanics (breaking surface gravity waves on a slope)
# Jakob number	Ja	\mathrm{Ja} = \frac{c_p (T_\mathrm{s} - T_\mathrm{sat}) }{\Delta H_{\mathrm{f}} }	chemistry (ratio of sensible to latent energy absorbed during liquid-vapor phase change)[18]
# Karlovitz number	Ka	\mathrm{Ka} = k t_c	turbulent combustion (characteristic flow time times flame stretch rate)
# Keulegan–Carpenter number	KC	\mathrm{K_C} = \frac{V\,T}{L}	fluid dynamics (ratio of drag force to inertia for a bluff object in oscillatory fluid flow)
# Knudsen number	Kn	\mathrm{Kn} = \frac {\lambda}{L}	gas dynamics (ratio of the molecular mean free path length to a representative physical length scale)
# Kt/V	Kt/V		medicine (hemodialysis and peritoneal dialysis treatment; dimensionless time)
# Kutateladze number	Ku	\mathrm{Ku} = \frac{U_h \rho_g^{1/2}}{\left({\sigma g (\rho_l - \rho_g)}\right)^{1/4}}	fluid mechanics (counter-current two-phase flow)[19]
# Laplace number	La	\mathrm{La} = \frac{\sigma \rho L}{\mu^2}	fluid dynamics (free convection within immiscible fluids; ratio of surface tension to momentum-transport)
# Lewis number	Le	\mathrm{Le} = \frac{\alpha}{D} = \frac{\mathrm{Sc}}{\mathrm{Pr}}	heat and mass transfer (ratio of thermal to mass diffusivity)
# Lift coefficient	CL	C_\mathrm{L} = \frac{L}{q\,S}	aerodynamics (lift available from an airfoil at a given angle of attack)
# Lockhart–Martinelli parameter	\chi	\chi = \frac{m_\ell}{m_g} \sqrt{\frac{\rho_g}{\rho_\ell}}	two-phase flow (flow of wet gases; liquid fraction)[20]
# Love numbers	h, k, l		geophysics (solidity of earth and other planets)
# Lundquist number	S	S = \frac{\mu_0LV_A}{\eta}	plasma physics (ratio of a resistive time to an Alfvén wave crossing time in a plasma)
# Mach number	M or Ma	 \mathrm{M} = \frac{{v}}{{v_\mathrm{sound}}}	gas dynamics (compressible flow; dimensionless velocity)
# Magnetic Reynolds number	Rm	\mathrm{R}_\mathrm{m} = \frac{U L}{\eta}	magnetohydrodynamics (ratio of magnetic advection to magnetic diffusion)
# Manning roughness coefficient	n		open channel flow (flow driven by gravity)[21]
# Marangoni number	Mg	\mathrm{Mg} = - {\frac{\mathrm{d}\sigma}{\mathrm{d}T}}\frac{L \Delta T}{\eta \alpha} 	fluid mechanics (Marangoni flow; thermal surface tension forces over viscous forces)
# Morton number	Mo	\mathrm{Mo} = \frac{g \mu_c^4 \, \Delta \rho}{\rho_c^2 \sigma^3}  	fluid dynamics (determination of bubble/drop shape)
# Nusselt number	Nu	\mathrm{Nu} =\frac{hd}{k}	heat transfer (forced convection; ratio of convective to conductive heat transfer)
# Ohnesorge number	Oh	 \mathrm{Oh} = \frac{ \mu}{ \sqrt{\rho \sigma L }} = \frac{\sqrt{\mathrm{We}}}{\mathrm{Re}} 	fluid dynamics (atomization of liquids, Marangoni flow)
# Péclet number	Pe	\mathrm{Pe} =  \frac{du\rho c_p}{k} = \mathrm{Re}\, \mathrm{Pr}	heat transfer (advection–diffusion problems; total momentum transfer to molecular heat transfer)
# Peel number	NP	N_\mathrm{P} = \frac{\text{Restoring force}}{\text{Adhesive force}}	coating (adhesion of microstructures with substrate)[22]
# Perveance	K	{K} = \frac{{I}}{{I_0}}\,\frac{{2}}{{\beta}^3{\gamma}^3} (1-\gamma^2f_e)	charged particle transport (measure of the strength of space charge in a charged particle beam)
# pH	\mathrm{pH}	\mathrm{pH} = - \log_{10}(a_{\textrm{H}^+}) = \log_{10}\left(\frac{1}{a_{\textrm{H}^+}}\right)	chemistry (the measure of the acidity or basicity of an aqueous solution)
# Pi	\pi	\pi = \frac{C}{d} \approx 3.14159	mathematics (ratio of a circle's circumference to its diameter)
# Pixel	px		digital imaging (smallest addressable unit)
# Poisson's ratio	\nu	\nu = -\frac{\mathrm{d}\varepsilon_\mathrm{trans}}{\mathrm{d}\varepsilon_\mathrm{axial}} 	elasticity (load in transverse and longitudinal direction)
# Porosity	\phi	\phi = \frac{V_\mathrm{V}}{V_\mathrm{T}}	geology, porous media (void fraction of the medium)
# Power factor	P/S		electronics (real power to apparent power)
# Power number	Np	 N_p = {P\over \rho n^3 d^5} 	electronics (power consumption by agitators; resistance force versus inertia force)
# Prandtl number	Pr	\mathrm{Pr} = \frac{\nu}{\alpha}  = \frac{c_p \mu}{k}	heat transfer (ratio of viscous diffusion rate over thermal diffusion rate)
# Prater number	β	\beta = \frac{-\Delta H_r D_{TA}^e C_{AS}}{\lambda^e T_s}	reaction engineering (ratio of heat evolution to heat conduction within a catalyst pellet)[23]
# Pressure coefficient	CP	C_p = {p - p_\infty \over \frac{1}{2} \rho_\infty V_\infty^2}	aerodynamics, hydrodynamics (pressure experienced at a point on an airfoil; dimensionless pressure variable)
# Q factor	Q	Q = 2 \pi f_r \frac{\text{Energy Stored}}{\text{Power Loss}}	physics, engineering (damping of oscillator or resonator; energy stored versus energy lost)
# Radian measure	rad	\text{arc length}/\text{radius}	mathematics (measurement of planar angles, 1 radian = 180/π degrees)
# Rayleigh number	Ra	\mathrm{Ra}_{x} = \frac{g \beta} {\nu \alpha} (T_s - T_\infin) x^3 	heat transfer (buoyancy versus viscous forces in free convection)
# Refractive index	n	n=\frac{c}{v}	electromagnetism, optics (speed of light in a vacuum over speed of light in a material)
# Relative density	RD	RD = \frac{\rho_\mathrm{substance}}{\rho_\mathrm{reference}}	hydrometers, material comparisons (ratio of density of a material to a reference material—usually water)
# Relative permeability	\mu_r	\mu_r = \frac{\mu}{\mu_0}	magnetostatics (ratio of the permeability of a specific medium to free space)
# Relative permittivity	\varepsilon_r	\varepsilon_{r} = \frac{C_{x}} {C_{0}}	electrostatics (ratio of capacitance of test capacitor with dielectric material versus vacuum)
# Reynolds number	Re	\mathrm{Re} = \frac{vL\rho}{\mu}	fluid mechanics (ratio of fluid inertial and viscous forces)[5]
# Richardson number	Ri	 \mathrm{Ri} = \frac{gh}{u^2} = \frac{1}{\mathrm{Fr}^2} 	fluid dynamics (effect of buoyancy on flow stability; ratio of potential over kinetic energy)[24]
# Rockwell scale	–		mechanical hardness (indentation hardness of a material)
# Rolling resistance coefficient	Crr	C_{rr} = \frac{F}{N_f} 	vehicle dynamics (ratio of force needed for motion of a wheel over the normal force)
# Roshko number	Ro	 \mathrm{Ro} = {f L^{2}\over \nu} =\mathrm{St}\,\mathrm{Re} 	fluid dynamics (oscillating flow, vortex shedding)
# Rossby number	Ro	\mathrm{Ro}=\frac{U}{Lf}	geophysics (ratio of inertial to Coriolis force)
# Rouse number	P or Z	\mathrm{P} = \frac{w_s}{\kappa u_*}	sediment transport (ratio of the sediment fall velocity and the upwards velocity of grain)

# Shape factor	H	H = \frac {\delta^*}{\theta}	boundary layer flow (ratio of displacement thickness to momentum thickness)
# Sherwood number	Sh	\mathrm{Sh} = \frac{K L}{D} 	mass transfer (forced convection; ratio of convective to diffusive mass transport)
# Shields parameter	\tau_* or \theta	\tau_{\ast} = \frac{\tau}{(\rho_s - \rho) g D}	sediment transport (threshold of sediment movement due to fluid motion; dimensionless shear stress)
# Sommerfeld number	S	 \mathrm{S} = \left( \frac{r}{c} \right)^2 \frac {\mu N}{P}	hydrodynamic lubrication (boundary lubrication)[26]
# Specific gravity	SG		(same as Relative density)
# Stanton number	St	\mathrm{St} = \frac{h}{c_p \rho V} = \frac{\mathrm{Nu}}{\mathrm{Re}\,\mathrm{Pr}} 	heat transfer and fluid dynamics (forced convection)
# Stefan number	Ste	\mathrm{Ste} = \frac{c_p \Delta T}{L}	phase change, thermodynamics (ratio of sensible heat to latent heat)
# Stokes number	Stk or Sk	\mathrm{Stk} = \frac{\tau U_o}{d_c}	particles suspensions (ratio of characteristic time of particle to time of flow)
# Strain	\epsilon	\epsilon = \cfrac{\partial{F}}{\partial{X}} - 1	materials science, elasticity (displacement between particles in the body relative to a reference length)
# Strouhal number	St or Sr	\mathrm{St} = {\omega L\over v} 	fluid dynamics (continuous and pulsating flow; nondimensional frequency)[27]
# Stuart number	N	 \mathrm{N} = \frac {B^2 L_{c} \sigma}{\rho U} = \frac{\mathrm{Ha}^2}{\mathrm{Re}} 	magnetohydrodynamics (ratio of electromagnetic to inertial forces)
# Taylor number	Ta	 \mathrm{Ta} = \frac{4\Omega^2 R^4}{\nu^2}	fluid dynamics (rotating fluid flows; inertial forces due to rotation of a fluid versus viscous forces)
# Transmittance	T	 T =  \frac{I}{I_0}	optics, spectroscopy (the ratio of the intensities of radiation exiting through and incident on a sample)
# Ursell number	U	\mathrm{U} = \frac{H\, \lambda^2}{h^3}	wave mechanics (nonlinearity of surface gravity waves on a shallow fluid layer)
# Vadasz number	Va	\mathrm{Va} = \frac{\phi\, \mathrm{Pr}}{\mathrm{Da}}	porous media (governs the effects of porosity \phi, the Prandtl number and the Darcy number on flow in a porous medium) [28]
# van 't Hoff factor	i	 i = 1 + \alpha (n - 1)	quantitative analysis (Kf and Kb)
# Wallis parameter	j*	j^* = R \left( \frac{\omega \rho}{\mu} \right)^\frac{1}{2}	multiphase flows (nondimensional superficial velocity)[29]
# Weaver flame speed number	Wea	\mathrm{Wea} = \frac{w}{w_\mathrm{H}} 100	combustion (laminar burning velocity relative to hydrogen gas)[30]

# Weissenberg number	Wi	\mathrm{Wi} = \dot{\gamma} \lambda 	viscoelastic flows (shear rate times the relaxation time)[31]
# Womersley number	\alpha	\alpha = R \left( \frac{\omega \rho}{\mu} \right)^\frac{1}{2}	biofluid mechanics (continuous and pulsating flows; ratio of pulsatile flow frequency to viscous effects)[32]

end # module

max_stars_repo_name: sabjohnso/DimensionalAnalysis

max_issues_repo_name: sabjohnso/DimensionalAnalysis

max_forks_repo_name: sabjohnso/DimensionalAnalysis


text:
using TkIntcode
using Test

@testset "TkIntcode.jl" begin

    # echo program
    let value = 123
        start(intcode("3,0,4,0,99", () -> value, x -> @test(x == value)))
    end

    # day5 part2 test1
    let code = "3,21,1008,21,8,20,1005,20,22,107,8,21,20,1006,20,31,1106,0,36,98,0,0,1002,21,125,20,4,20,1105,1,46,104,999,1105,1,46,1101,1000,1,20,4,20,1105,1,46,98,99"
        for value in 1:7
            start(intcode(code, () -> value, x -> @test(x == 999)))
        end
        start(intcode(code, () -> 8, x -> @test(x == 1000)))
        start(intcode(code, () -> 9, x -> @test(x == 1001)))
        start(intcode(code, () -> 10, x -> @test(x == 1001)))
    end

end

ext: jl


text:
# QuEST.jl/test/runtests.jl
#
# Authors:
#  - Dirk Oliver Theis, Ketita Labs
#
# Copyright (c) 2020 Ketita Labs oü, Tartu, Estonia
#
# MIT License
#
# Copyright (c) 2020 Ketita Labs oü, Tartu, Estonia
#
# Permission is hereby granted, free of charge, to any person
# obtaining a copy of this software and associated documentation files
# (the "Software"), to deal in the Software without restriction,
# including without limitation the rights to use, copy, modify, merge,
# publish, distribute, sublicense, and/or sell copies of the Software,
# and to permit persons to whom the Software is furnished to do so,
# subject to the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
# BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

using QuEST
using Test

function test_env_and_qureg()

    function test_2quregs(q ::QuEST.Qureg, env ::QuEST.QuESTEnv)
        q1 = QuEST.createCloneQureg(q,env)
        q2 = QuEST.createCloneQureg(q,env)
        @warn "need more here"
    end
    function test_qureg(env ::QuEST.QuESTEnv)
        qureg = QuEST.createQureg(8,env)
        test_2quregs(qureg,env)
        QuEST.destroyQureg(qureg,env)
        @warn "need more here"
    end

    function test_densiQureg(env ::QuEST.QuESTEnv)
        qureg = QuEST.createDensityQureg(8,env)
        test_2quregs(qureg,env)
        QuEST.destroyQureg(qureg,env)
        @warn "need more here"
    end

    function test_matrices()
        M = QuEST.createComplexMatrixN(3)
        QuEST.fill_ComplexMatrix!(M, (k,ℓ) -> exp(im*(k+ℓ))/sqrt(8) )
        QuEST.destroyComplexMatrixN(M)
        @warn "need more here"
    end

    env = QuEST.createQuESTEnv()
    test_qureg(env)
    test_densiQureg(env)
    QuEST.destroyQuESTEnv(env)

    test_matrices()
end

#
@testset "QuEST.jl" begin
    test_env_and_qureg()
end

# EOF


text:
# Container for a DataTable row
immutable DataTableRow{T <: AbstractDataTable}
    dt::T
    row::Int
end

function Base.getindex(r::DataTableRow, idx::AbstractArray)
    return DataTableRow(r.dt[idx], r.row)
end

function Base.getindex(r::DataTableRow, idx::Any)
    return r.dt[r.row, idx]
end

function Base.setindex!(r::DataTableRow, value::Any, idx::Any)
    return setindex!(r.dt, value, r.row, idx)
end

Base.names(r::DataTableRow) = names(r.dt)
_names(r::DataTableRow) = _names(r.dt)

Base.view(r::DataTableRow, c) = DataTableRow(r.dt[[c]], r.row)

index(r::DataTableRow) = index(r.dt)

Base.length(r::DataTableRow) = size(r.dt, 2)

Base.endof(r::DataTableRow) = size(r.dt, 2)

Base.collect(r::DataTableRow) = Tuple{Symbol, Any}[x for x in r]

Base.start(r::DataTableRow) = 1

Base.next(r::DataTableRow, s) = ((_names(r)[s], r[s]), s + 1)

Base.done(r::DataTableRow, s) = s > length(r)

Base.convert(::Type{Array}, r::DataTableRow) = convert(Array, r.dt[r.row,:])

# hash column element
Base.@propagate_inbounds hash_colel(v::AbstractArray, i, h::UInt = zero(UInt)) = hash(v[i], h)
Base.@propagate_inbounds hash_colel{T<:Nullable}(v::AbstractArray{T}, i, h::UInt = zero(UInt)) =
    isnull(v[i]) ? h + Base.nullablehash_seed : hash(unsafe_get(v[i]), h)
Base.@propagate_inbounds hash_colel{T}(v::NullableArray{T}, i, h::UInt = zero(UInt)) =
    isnull(v, i) ? h + Base.nullablehash_seed : hash(v.values[i], h)
Base.@propagate_inbounds hash_colel{T}(v::AbstractCategoricalArray{T}, i, h::UInt = zero(UInt)) =
    hash(CategoricalArrays.index(v.pool)[v.refs[i]], h)
Base.@propagate_inbounds function hash_colel{T}(v::AbstractNullableCategoricalArray{T}, i, h::UInt = zero(UInt))
    ref = v.refs[i]
    ref == 0 ? h + Base.nullablehash_seed : hash(CategoricalArrays.index(v.pool)[ref], h)
end

# hash of DataTable rows based on its values
# so that duplicate rows would have the same hash
# table columns are passed as a tuple of vectors to ensure type specialization
rowhash(cols::Tuple{AbstractVector}, r::Int, h::UInt = zero(UInt))::UInt =
    hash_colel(cols[1], r, h)
function rowhash(cols::Tuple{Vararg{AbstractVector}}, r::Int, h::UInt = zero(UInt))::UInt
    h = hash_colel(cols[1], r, h)
    rowhash(Base.tail(cols), r, h)
end

Base.hash(r::DataTableRow, h::UInt = zero(UInt)) =
    rowhash(ntuple(i -> r.dt[i], ncol(r.dt)), r.row, h)

# comparison of DataTable rows
# only the rows of the same DataTable could be compared
# rows are equal if they have the same values (while the row indices could differ)
# returns Nullable{Bool}
# if all non-null values are equal, but there are nulls, returns null
@compat(Base.:(==))(r1::DataTableRow, r2::DataTableRow) = isequal(r1, r2)

function Base.isequal(r1::DataTableRow, r2::DataTableRow)
    isequal_row(r1.dt, r1.row, r2.dt, r2.row)
end

# internal method for comparing the elements of the same data table column
isequal_colel(col::AbstractArray, r1::Int, r2::Int) =
    (r1 == r2) || isequal(Base.unsafe_getindex(col, r1), Base.unsafe_getindex(col, r2))

isequal_colel(a::Any, b::Any) = isequal(a, b)
isequal_colel(a::Nullable, b::Any) = !isnull(a) & isequal(unsafe_get(a), b)
isequal_colel(a::Any, b::Nullable) = isequal_colel(b, a)
isequal_colel(a::Nullable, b::Nullable) = isequal(a, b)

# table columns are passed as a tuple of vectors to ensure type specialization
isequal_row(cols::Tuple{AbstractVector}, r1::Int, r2::Int) =
    isequal_colel(cols[1][r1], cols[1][r2])
isequal_row(cols::Tuple{Vararg{AbstractVector}}, r1::Int, r2::Int) =
    isequal_colel(cols[1][r1], cols[1][r2]) && isequal_row(Base.tail(cols), r1, r2)

isequal_row(cols1::Tuple{AbstractVector}, r1::Int, cols2::Tuple{AbstractVector}, r2::Int) =
    isequal_colel(cols1[1][r1], cols2[1][r2])
isequal_row(cols1::Tuple{Vararg{AbstractVector}}, r1::Int,
            cols2::Tuple{Vararg{AbstractVector}}, r2::Int) =
    isequal_colel(cols1[1][r1], cols2[1][r2]) &&
        isequal_row(Base.tail(cols1), r1, Base.tail(cols2), r2)

function isequal_row(dt1::AbstractDataTable, r1::Int, dt2::AbstractDataTable, r2::Int)
    if dt1 === dt2
        if r1 == r2
            return true
        end
    elseif !(ncol(dt1) == ncol(dt2))
        throw(ArgumentError("Rows of the tables that have different number of columns cannot be compared. Got $(ncol(dt1)) and $(ncol(dt2)) columns"))
    end
    @inbounds for (col1, col2) in zip(columns(dt1), columns(dt2))
        isequal_colel(col1[r1], col2[r2]) || return false
    end
    return true
end

# lexicographic ordering on DataTable rows, null > !null
function Base.isless(r1::DataTableRow, r2::DataTableRow)
    (ncol(r1.dt) == ncol(r2.dt)) ||
        throw(ArgumentError("Rows of the data tables that have different number of columns cannot be compared ($(ncol(dt1)) and $(ncol(dt2)))"))
    @inbounds for i in 1:ncol(r1.dt)
        x = r1.dt[i][r1.row]
        y = r2.dt[i][r2.row]
        isnullx = _isnull(x)
        isnully = _isnull(y)
        (isnullx != isnully) && return isnully # null > !null
        if !isnullx
            v1 = unsafe_get(x)
            v2 = unsafe_get(y)
            isless(v1, v2) && return true
            !isequal(v1, v2) && return false
        end
    end
    return false
end


text:
@testset "volume-volume" begin

    ### Box properties
    box_rad = 0.05
    box_density = 400.0
    w_z_0 = 1.14
    i_prop_compliant = InertiaProperties(box_density)
    c_prop_compliant = ContactProperties(Ē=1.0e6)
    eM_box_compliant = as_tet_eMesh(eMesh_box(box_rad))

    mech_scen = MechanismScenario()  # n_quad_rule=2)

    ### Add planes and boxes
    nt_plane  = add_contact!(     mech_scen, "plane", as_tet_eMesh(eMesh_half_plane()),   c_prop=c_prop_compliant)
    nt_body_1 = add_body_contact!(mech_scen, "box_1", eM_box_compliant, i_prop=i_prop_compliant, c_prop=c_prop_compliant)

    ### Friction
    add_friction_regularize!(mech_scen, nt_plane.id,  nt_body_1.id, μd=0.0, χ=0.0, n_quad_rule=2)

    finalize!(mech_scen)
    set_state_spq!(mech_scen, nt_body_1.joint, trans=SVector(0.0, 0.0, 2*box_rad), w=SVector(0.0, 0.0, w_z_0))

    ### Run forward dynamics
    t_final = 5.0e-0
    rr = Radau_for_MechanismScenario(mech_scen)
    rr.step.h_max = 0.05
    data_time, data_state = integrate_scenario_radau(rr, t_final=t_final)

    @test data_state[end, 9] ≈ w_z_0
end


text:
model=GeoPhyInv.Gallery.Seismic(:acou_homo1)
GeoPhyInv.Models.Seismic_addon!(model, randn_perc=0.1, fields=[:χvp,:χρ])
acqgeom =GeoPhyInv.Acquisition.Geom_circ(nss=4,nr=100,rad=[990.,990.]);
acqgeom =GeoPhyInv.Acquisition.Geom_circ(nss=4,nr=100,rad=[0.,200.]);
acqsrc=GeoPhyInv.Acquisition.Src_fixed_mod(acqgeom.nss,1,[:P],mod=model, nλ=3, tmaxfrac=0.4)

for sflags in [[1,-1],[2,-2]]
	pa=GeoPhyInv.Fdtd.Param(born_flag=false,npw=1, tgridmod=acqsrc.tgrid,
	#	abs_trbl=[:null],
		gmodel_flag=false,
		sflags=[sflags[1]],
		snaps_flag=true,
		verbose=true,
		backprop_flag=1,
		illum_flag=true,acqgeom=[acqgeom], acqsrc=[acqsrc],
		model=model);

	GeoPhyInv.Fdtd.mod!(pa);
	rec1=deepcopy(pa.c.data[1])

	# change source flag and update wavelets in pa
	pa.c.sflags=[sflags[2]];
	GeoPhyInv.Fdtd.update_acqsrc!(pa,[acqsrc])
	pa.c.backprop_flag=-1 # do backpropagation

	GeoPhyInv.Fdtd.mod!(pa);
	rec2=deepcopy(pa.c.data[1])

	# time reverse
	GeoPhyInv.Data.TD_tr!(rec2);

	# compare results
	# least-squares misfit
	paerr=GeoPhyInv.Data.P_misfit(rec1, rec2)
	err = GeoPhyInv.Data.func_grad!(paerr)

	# normalized error
	error = err[1]./paerr.ynorm

	# desired accuracy?
	@test error<1e-20
end


text:
using Revise
using IsopycnalSurfaces, Test

pz = collect(0.:500.:4000.) # pressure levels
nz = length(pz)
θz = collect(range(20,stop=10,length=nz))
Sz = collect(range(36,stop=35,length=nz))

println("Using the EOS from MITgcmTools.jl:")
display(sigma2column(θz,Sz,pz,"EOS94")) # using the EOS from MITgcmTools.jl
println("")
println("Default, using the EOS from PhysOcean.jl/EOS80.jl:")
display(sigma2column(θz,Sz,pz))  # default, using the EOS from PhysOcean.jl/EOS80.jl


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

"""
    month_cnv(number[, shor=true, up=true, low=true]) -> month_name
    month_cnv(name) -> number

### Purpose ###

Convert between a month English name and  the equivalent number.

### Explanation ###

For example, converts from "January" to 1  or vice-versa.

### Arguments ###

The functions has two methods, one with numeric input (and three possible
boolean keywords) and the other one with string input.

Numeric input arguments:

* `number`: the number of the month to be converted to month name.
* `short` (optional boolean keyword): if true, the abbreviated (3-character)
  name of the month will be returned, e.g. "Apr" or "Oct".  Default is false.
* `up` (optional boolean keyword): if true, the name of the month will be all in
  upper case, e.g. "APRIL" or "OCTOBER".  Default is false.
* `low` (optional boolean keyword): if true, the name of the month will be all
  in lower case, e.g. "april" or "october".  Default is false.

String input argument:

* `name`: month name to be converted to month number.

### Output ###

The month name or month number, depending on the input.  For numeric input, the
format of the month name is influenced by the optional keywords.

### Example ###

```jldoctest
julia> using AstroLib

julia> month_cnv.(["janua", "SEP", "aUgUsT"])
3-element Array{Int64,1}:
 1
 9
 8

julia> month_cnv.([2, 12, 6], short=true, low=true)
3-element Array{String,1}:
 "feb"
 "dec"
 "jun"
```

"""
function month_cnv(number::Integer; short::Bool=false,
                   up::Bool=false, low::Bool=false)
    if short
        name = Dates.ENGLISH.months_abbr[number]
    else
        name = Dates.ENGLISH.months[number]
    end
    if up
        name = uppercase(name)
    elseif low
        name = lowercase(name)
    end
    return name
end

function month_cnv(name::AbstractString)
    name = strip(name)
    if length(name) >= 3
        output::Integer = getfield(Dates, Symbol(uppercasefirst(lowercase(name[1:3]))))
        return output
    else
        # Do the same as original AstroLib: return -1 for input string shorter
        # than 3.
        return -1
    end
end


text:
#=
# A struct containing configuration information
# This is the most general Finch info that should be (almost)problem indepenent.
# Premade configs can be loaded for different classes of problems.
# Values can be set through various macros or automatically by problem specification.
=#
if !@isdefined(IRREGULAR)
    include("finch_constants.jl");
end

mutable struct Finch_config
    # Domain
    dimension::Int          # 1,2,3
    geometry::String        # square, irregular
    mesh_type::String       # unstructured, tree, uniform grid

    # FEM details
    solver_type::String     # cg, dg, fv
    trial_function::String  # Legendre
    test_function::String   # same as above
    elemental_nodes::String # uniform, gauss, lobatto (higher order node distribution within elements)
    quadrature::String      # uniform, gauss, lobatto (similar to above)
    p_adaptive::Bool        # Do adaptive p-refinement?
    basis_order_min::Int    # minimum order to use in p-refinement, or if p_adaptive is false
    basis_order_max::Int    # maximum order

    # Other solver details
    linear::Bool            # Is the equation linear?
    t_adaptive::Bool        # Do adaptive t_refinement?
    stepper::String         # Euler-explicit/implicit, RK4, LSRK4, etc. Type of time stepper to use
    linalg_matrixfree::Bool # Use matrix free methods?
    linalg_matfree_max::Int # max iters for matrix free
    linalg_matfree_tol::Float64 # tolerance for matrix free
    linalg_backend::String  # default, petsc, ?? (What to use for linear algebra)

    # Output
    output_format::String   # VTK, raw, custom (format for storing solutions)

    # Parallel details
    use_mpi::Bool           # Is MPI available?
    num_procs::Int;         # number of processes
    proc_rank::Int;         # this proccess rank
    num_threads::Int;       # number of available threads
    num_partitions::Int;    # number of mesh partitions
    partition_index::Int;   # this process's partition

    # Constructor builds a default config.
    Finch_config() = new(
        1,
        SQUARE,
        UNIFORM_GRID,
        CG,
        LEGENDRE,
        LEGENDRE,
        LOBATTO,
        GAUSS,
        false,
        1,
        1,
        true,
        false,
        EULER_IMPLICIT,
        false,
        1,
        1.0,
        DEFAULT_SOLVER,
        VTK,
        false,
        1,
        0,
        1,
        1,
        0
    );
end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
@field ColorBuffer
@field DepthBuffer

@field Model = eye(Mat4f0)
@field Renderlist = []
@field NativeWindow = nothing
"""
Holds general attribute of the scene/window, like camera uniforms and lights.
"""
@field Scene = nothing
@field Canvas = nothing

abstract type AbstractScene <: ReactiveComposable end
abstract type AbstractCanvas <: ReactiveComposable end

@reactivecomposed type JLCanvas <: AbstractCanvas
    Area
    DepthBuffer
    ColorBuffer
end

function default(::Type{DepthBuffer}, canvas::Partial{<: AbstractCanvas})
    w, h = widths(get(canvas, Area))
    ones(Float32, w, h)
end

function default(::Type{ColorBuffer}, canvas::Partial{<: AbstractCanvas})
    w, h = widths(get(canvas, Area))
    (zeros(RGBA{Float32}, w, h), )
end

@field Ambient = Vec3f0(0.1)
@field Diffuse = Vec3f0(0.0)
@field DiffusePower = 1f0
@field Specular = Vec3f0(0.3)
@field SpecularPower = 10f0

@composed type Light
    Position::Vec3f0
    Ambient::Vec3f0
    Diffuse::Vec3f0
    DiffusePower::Float32
    Specular::Vec3f0
    SpecularPower::Float32
end

@composed type Window <: AbstractWindow
    <: WindowEvents
    NativeWindow
    Camera
    Renderlist
    Scene
    Color
    Light # TODO Move into Scene, but can't since struct of structs isn't supported very well yet
end
function default(::Type{Scene}, p::Partial{<: AbstractWindow})
    cam = get!(p, Camera)
    scene = SceneUniforms(cam)
    for field in FieldTraits.Fields(scene)
        if haskey(cam, field)
            FieldTraits.link!(field, cam => scene)
        end
    end
    return scene
end

Base.isopen(window::AbstractWindow) = window[Open]

function getmultiple!(field, parent, tail...)
    i = start(tail)
    head = parent
    while true
        haskey(head, field) && return get!(head, field)
        if done(tail, i)
            def = default(field, parent)
            parent[field] = def
            return def
        end
        head, i = next(tail, i)
    end
end

# function default{AC <: AbstractScene}(::Type{AC}, p::Partial{Window})
#     cam = get!(p, Camera)
#     AC(map(field-> getmultiple!(field, p, cam), Fields(AC)))
# end

# Basic functions a window should defined
function show!(window::AbstractWindow)
    error("Not implemented for $(typeof(window))")
end
function destroy!(window::AbstractWindow)
    error("Not implemented for $(typeof(window))")
end
function swapbuffers!(window::AbstractWindow)
    error("Not implemented for $(typeof(window))")
end
function renderloop(window::AbstractWindow)
    error("Not implemented for $(typeof(window))")
end


text:
using ITensors
using ITensorTDVP
using Observers

function heisenberg(N)
  os = OpSum()
  for j in 1:(N - 1)
    os += 0.5, "S+", j, "S-", j + 1
    os += 0.5, "S-", j, "S+", j + 1
    os += "Sz", j, "Sz", j + 1
  end
  return os
end

N = 10
cutoff = 1e-12
tau = 0.1
ttotal = 1.0

s = siteinds("S=1/2", N; conserve_qns=true)
H = MPO(heisenberg(N), s)

function step(; sweep, bond, half_sweep)
  if bond == 1 && half_sweep == 2
    return sweep
  end
  return nothing
end

function current_time(; current_time, bond, half_sweep)
  if bond == 1 && half_sweep == 2
    return current_time
  end
  return nothing
end

function measure_sz(; psi, bond, half_sweep)
  if bond == 1 && half_sweep == 2
    return expect(psi, "Sz"; sites=N ÷ 2)
  end
  return nothing
end

function return_state(; psi, bond, half_sweep)
  if bond == 1 && half_sweep == 2
    return psi
  end
  return nothing
end

obs = Observer(
  "steps" => step, "times" => current_time, "psis" => return_state, "Sz" => measure_sz
)

psi = MPS(s, n -> isodd(n) ? "Up" : "Dn")
psi_f = tdvp(
  H,
  -im * ttotal,
  psi;
  time_step=-im * tau,
  cutoff,
  outputlevel=1,
  normalize=false,
  (observer!)=obs,
)

res = results(obs)
steps = res["steps"]
times = res["times"]
psis = res["psis"]
Sz = res["Sz"]

println("\nResults")
println("=======")
for n in 1:length(steps)
  print("step = ", steps[n])
  print(", time = ", round(times[n]; digits=3))
  print(", |⟨ψⁿ|ψⁱ⟩| = ", round(abs(inner(psis[n], psi)); digits=3))
  print(", |⟨ψⁿ|ψᶠ⟩| = ", round(abs(inner(psis[n], psi_f)); digits=3))
  print(", ⟨Sᶻ⟩ = ", round(Sz[n]; digits=3))
  println()
end


text:
using Test

@testset "InternalForces" begin

    include(string(@__DIR__) * "/InternalForces/" * "InternalForcesTest1.jl")
    include(string(@__DIR__) * "/InternalForces/" * "InternalForcesTest2.jl")
    include(string(@__DIR__) * "/InternalForces/" * "InternalForcesTest3.jl")
    include(string(@__DIR__) * "/InternalForces/" * "InternalForcesTest4.jl")
    include(string(@__DIR__) * "/InternalForces/" * "InternalForcesTest5.jl")

end


text:
export vec_data

abstract type AbstractDerivative end

struct RealDerivative{NT,V} <: AbstractDerivative
    fields::NT
    vectorised_data::V
end

@inline Base.propertynames(s::RealDerivative) = propertynames(getfield(s, :fields))
@inline Base.getindex(s::RealDerivative, val::Symbol) =
    getproperty(s, val)
@inline function Base.getproperty(s::RealDerivative, val::Symbol)
    val===:tuple_all_weights && return vec_data(s)
    return getproperty(getfield(s, :fields), val)
end
@inline vec_data(s::RealDerivative) = getfield(s, :vectorised_data)

struct WirtingerDerivative{R,V} <: AbstractDerivative
    r_derivatives::R
    c_derivatives::R
    vectorised_data::V
end

vec_data(s::WirtingerDerivative)  = s.vectorised_data
Base.real(s::WirtingerDerivative) = s.r_derivatives
Base.imag(s::WirtingerDerivative) = s.c_derivatives

@inline function Base.getproperty(s::WirtingerDerivative, val::Symbol)
    val===:tuple_all_weights && return vec_data(s)
    return getfield(s, val)
end

function RealDerivative(net::NeuralNetwork)
    vec    = Vector{out_type(net)}()
    i, fields = weight_tuple(net, fieldnames(typeof(net)), vec)
    return RealDerivative(fields, [vec])
end

function WirtingerDerivative(net::NeuralNetwork)
    vec         = Vector{out_type(net)}()
    i, fields_r = weight_tuple(net, fieldnames(typeof(net)), vec)
    i, fields_c = weight_tuple(net, fieldnames(typeof(net)), vec, i+1)
    return WirtingerDerivative(fields_r, fields_c, [vec])
end

Base.show(io::IO, der::RealDerivative) = begin
    pn = propertynames(der)
    str = "{"
    for fn=pn[1:end-1]
        str *= ":$fn, "
    end
    str *= ":" * string(last(pn))*" }"
    print(io,
    "RealDerivative with fields: ", str)
end

Base.show(io::IO, ::MIME"text/plain", der::RealDerivative) = begin
    pn = propertynames(der)
    str = "{"
    for fn=pn[1:end-1]
        str *= ":$fn, "
    end
    str *= ":" * string(last(pn))*" }"
    print(io,
    "RealDerivative with fields: ", str)
end

Base.show(io::IO, der::WirtingerDerivative) = begin
    pn = propertynames(real(der))
    str = "{"
    for fn=pn[1:end-1]
        str *= ":$fn, "
    end
    str *= ":" * string(last(pn))*" }"
    print(io,
    "WirtingerDerivative with fields: ", str)
end

Base.show(io::IO, ::MIME"text/plain", der::WirtingerDerivative) = begin
    pn = propertynames(real(der))
    str = "{"
    for fn=pn[1:end-1]
        str *= ":$fn, "
    end
    str *= ":" * string(last(pn))*" }"
    print(io,
    "WirtingerDerivative with fields: ", str)
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    IsotropicBall(radius, [metric])

A ball with single `radius` and `metric`.
Default metric is `Euclidean()`.
"""
struct IsotropicBall{T,M} <: MetricBall
  radius::T
  metric::M
end

function IsotropicBall(radius)
  metric = Euclidean()
  T = typeof(radius)
  M = typeof(metric)
  IsotropicBall{T,M}(radius, metric)
end

metric(ball::IsotropicBall) = ball.metric

radii(ball::IsotropicBall)  = SVector(ball.radius,)

function Base.show(io::IO, ball::IsotropicBall{T}) where {T}
  r, d = ball.radius, ball.metric
  print(io, "IsotropicBall{$T}($r, $d)")
end


text:
"""
The 2-W transformer model uses an equivalent circuit assuming the impedance is on the High Voltage Side of the transformer.
The model allocates the iron losses and magnetezing suceptance to the primary side
"""
struct Transformer2W <: Branch
    name::String
    available::Bool
    connectionpoints::From_To_Bus
    r::Float64 #[pu]
    x::Float64 #[pu]
    primaryshunt::Float64 #[pu]
    rate::Union{Nothing,Float64} #[MVA]
end

Transformer2W(; name = "init",
                available = false,
                connectionpoints = (from = Bus(), to =Bus()),
                r = 0.0,
                x = 0.0,
                primaryshunt = 0.0,
                rate = nothing
            ) = Transformer2W(name, available, connectionpoints, r, x, primaryshunt, rate)

struct TapTransformer <: Branch
    name::String
    available::Bool
    connectionpoints::From_To_Bus
    r::Float64 #[pu]
    x::Float64 #[pu]
    primaryshunt::Float64 #[pu]
    tap::Float64 # [0 - 2]
    rate::Union{Float64,Nothing} #[MVA]
end

TapTransformer(; name = "init",
                available = false,
                connectionpoints = (from=Bus(), to=Bus()),
                r = 0.0,
                x = 0.0,
                primaryshunt = 0.0,
                tap = 1.0,
                rate = nothing
            ) = TapTransformer(name, available, connectionpoints, r, x, primaryshunt, tap, rate)

#=
struct Transformer3W <: Branch
    name::String
    available::Bool
    transformer::Transformer2W
    line::Line
end

Transformer3W(; name = "init",
                available = false,
                transformer = Transformer2W(),
                line = Line()
            ) = Transformer3W(name, available, transformer, line)
=#

struct PhaseShiftingTransformer <: Branch
    name::String
    available::Bool
    connectionpoints::From_To_Bus
    r::Float64 #[pu]
    x::Float64 #[pu]
    primaryshunt::Float64 #[pu]
    tap::Float64 #[0 - 2]
    α::Float64 # [radians]
    rate::Float64 #[MVA]
end

PhaseShiftingTransformer(; name = "init",
                available = false,
                connectionpoints = (from=Bus(), to=Bus()),
                r = 0.0,
                x = 0.0,
                primaryshunt=0.0,
                tap = 1.0,
                α = 0.0,
                rate = 0.0
            ) = PhaseShiftingTransformer(name, available, connectionpoints, r, x, primaryshunt, tap, α, rate)


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

include("model_AccessUri.jl")
include("model_ApiError.jl")
include("model_ApiErrorBase.jl")
include("model_CloudError.jl")
include("model_CreationData.jl")
include("model_Disk.jl")
include("model_DiskAccess.jl")
include("model_DiskAccessList.jl")
include("model_DiskAccessProperties.jl")
include("model_DiskAccessUpdate.jl")
include("model_DiskEncryptionSet.jl")
include("model_DiskEncryptionSetList.jl")
include("model_DiskEncryptionSetUpdate.jl")
include("model_DiskEncryptionSetUpdateProperties.jl")
include("model_DiskList.jl")
include("model_DiskProperties.jl")
include("model_DiskSku.jl")
include("model_DiskUpdate.jl")
include("model_DiskUpdateProperties.jl")
include("model_Encryption.jl")
include("model_EncryptionSetIdentity.jl")
include("model_EncryptionSetProperties.jl")
include("model_EncryptionSettingsCollection.jl")
include("model_EncryptionSettingsElement.jl")
include("model_EncryptionType.jl")
include("model_GrantAccessData.jl")
include("model_ImageDiskReference.jl")
include("model_InnerError.jl")
include("model_KeyVaultAndKeyReference.jl")
include("model_KeyVaultAndSecretReference.jl")
include("model_NetworkAccessPolicy.jl")
include("model_PrivateEndpoint.jl")
include("model_PrivateEndpointConnection.jl")
include("model_PrivateEndpointConnectionProperties.jl")
include("model_PrivateEndpointConnectionProvisioningState.jl")
include("model_PrivateEndpointServiceConnectionStatus.jl")
include("model_PrivateLinkResource.jl")
include("model_PrivateLinkResourceListResult.jl")
include("model_PrivateLinkResourceProperties.jl")
include("model_PrivateLinkServiceConnectionState.jl")
include("model_Resource.jl")
include("model_ShareInfoElement.jl")
include("model_Snapshot.jl")
include("model_SnapshotList.jl")
include("model_SnapshotProperties.jl")
include("model_SnapshotSku.jl")
include("model_SnapshotUpdate.jl")
include("model_SnapshotUpdateProperties.jl")
include("model_SourceVault.jl")


text:
using CSV
using JuMP
using CPLEX
using DataFrames
using Random

include("functions.jl")

dataSet = "kidney"
dataFolder = "../data/"
resultsFolder = "../res/"

# Create the features tables (or load them if they already exist)
# Note: each line corresponds to an individual, the 1st column of each table contain the class
# Details:
# - read the file ./data/kidney.csv
# - save the features in ./data/kidney_test.csv and ./data/kidney_train.csv
train, test = createFeatures(dataFolder, dataSet)

# Create the rules (or load them if they already exist)
# Note: each line corresponds to a rule, the first column corresponds to the class
# Details:
# - read the file ./data/kidney_train.csv
# - save the rules in ./res/kidney_rules.csv
rules = createRules(dataSet, resultsFolder, train)

# Order the rules (limit the resolution to 300 seconds)
# Details:
# - read the file ./data/kidney_rules.csv
# - save the rules in ./res/kidney_ordered_rules.csv
timeLimitInSeconds = 300
orderedRules = sortRules(dataSet, resultsFolder, train, rules, timeLimitInSeconds)

println("-- Train results")
showStatistics(orderedRules, train)

println("-- Test results")
showStatistics(orderedRules, test)


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    atol(T)

Absolute tolerance used in algorithms for approximate
comparison with numbers of type `T`. It is used in the
source code in calls to the [`isapprox`](@ref) function:

```julia
isapprox(a::T, b::T, atol=tol(T))
```
"""
atol(::Type{Float64}) = 1e-10
atol(::Type{Float32}) = 1f-6


text:
"""
Handles the non-relativistic NBody problem for planets,
and other celestial bodies.
"""
module NBody

using Reexport
@reexport using ..AstrodynamicsCore

include("../Misc/DocStringExtensions.jl")
include("../Misc/UnitfulAliases.jl")

using StaticArrays: SVector, @SVector, SMatrix, @SMatrix

export Body, NBodySystem, system_energy,
       system_angular_momentum, promote, convert,
       Float16, Float32, Float64, BigFloat,
       length, getindex

include("NBodyStates.jl")
include("NBodyCalculations.jl")

end


text:
mutable struct WordVectors{S<:AbstractString, T<:Real, H<:Integer}
    vocab::Vector{S} # vocabulary
    vectors::Array{T, 2} # the vectors computed from GloVe
    vocab_hash::Dict{S, H}
end

function WordVectors(vocab::AbstractArray{S,1},
                     vectors::AbstractArray{T,2}) where {S<: AbstractString, T<:Real}
    length(vocab) == size(vectors, 2) ||
        throw(DimensionMismatch("Dimension of vocab and vectors are inconsistent."))
    vocab_hash = Dict{S, Int}()
    for (i, word) in enumerate(vocab)
        vocab_hash[word] = i
    end
    WordVectors(vocab, vectors, vocab_hash)
end

function Base.show(io::IO, wv::WordVectors{S,T}) where {S,T}
    len_vecs, num_words = size(wv.vectors)
    print(io, "WordVectors $(num_words) words, $(len_vecs)-element $(T) vectors")
end

"""
    vocabulary(wv)

Return the vocabulary as a vector of words of the WordVectors `wv`.
"""
vocabulary(wv::WordVectors) = wv.vocab

"""
    in_vocabulary(wv, word)

Return `true` if `word` is part of the vocabulary of the WordVector `wv` and
`false` otherwise.
"""
in_vocabulary(wv::WordVectors, word::AbstractString) = word in wv.vocab

"""
    size(wv)

Return the word vector length and the number of words as a tuple.
"""
size(wv::WordVectors) = size(wv.vectors)

"""
    index(wv, word)

Return the index of `word` from the WordVectors `wv`.
"""
index(wv::WordVectors, word) = wv.vocab_hash[word]

"""
    get_vector(wv, word [; oov=false, oov_key="<unk>")

Returns the vector representation of `word` from the WordVectors `wv`.
If `oov` is `true`, the vector corresponding to the key `oov_key` is
returned for out-of-vocabulary words.
"""
get_vector(wv::WordVectors,
           word;
           oov::Bool=false,
           oov_key::String="<unk>") = begin
    if !oov
        idx = wv.vocab_hash[word];
        return wv.vectors[:, idx]
    else
        idx = get(wv.vocab_hash, word, wv.vocab_hash[oov_key])
        return wv.vectors[:, idx]
    end
end

"""
    cosine(wv, word, n=10)

Return the position of `n` (by default `n = 10`) neighbors of `word` and their
cosine similarities.
"""
function cosine(wv::WordVectors, word, n=10)
    metrics = wv.vectors'*get_vector(wv, word)
    topn_positions = sortperm(metrics[:], rev = true)[1:n]
    topn_metrics = metrics[topn_positions]
    return topn_positions, topn_metrics
end

"""
    similarity(wv, word1, word2)

Return the cosine similarity value between two words `word1` and `word2`.
"""
function similarity(wv::WordVectors, word1, word2)
    return get_vector(wv, word1)'*get_vector(wv, word2)
end

"""
    cosine_similar_words(wv, word, n=10)

Return the top `n` (by default `n = 10`) most similar words to `word`
from the WordVectors `wv`.
"""
function cosine_similar_words(wv::WordVectors, word, n=10)
    indx, metr = cosine(wv, word, n)
    return vocabulary(wv)[indx]
end

"""
    analogy(wv, pos, neg, n=5)

Compute the analogy similarity between two lists of words. The positions
and the similarity values of the top `n` similar words will be returned.
For example,
`king - man + woman = queen` will be
`pos=[\"king\", \"woman\"], neg=[\"man\"]`.
"""
function analogy(wv::WordVectors{S,T,H}, pos::AbstractArray, neg::AbstractArray, n= 5
                ) where {S<:AbstractString, T<:Real, H<:Integer}
    m, n_vocab = size(wv)
    n_pos = length(pos)
    n_neg = length(neg)
    anal_vecs = Matrix{T}(undef, m, n_pos + n_neg)

    for (i, word) in enumerate(pos)
        anal_vecs[:,i] = get_vector(wv, word)
    end
    for (i, word) in enumerate(neg)
        anal_vecs[:,i+n_pos] = -get_vector(wv, word)
    end
    mean_vec = mean(anal_vecs, dims=2)
    metrics = wv.vectors'*mean_vec
    top_positions = sortperm(metrics[:], rev = true)[1:n+n_pos+n_neg]
    for word in [pos;neg]
        idx = index(wv, word)
        loc = findfirst(x->x==idx, top_positions)
        if loc != nothing
            splice!(top_positions, loc)
        end
    end
    topn_positions = top_positions[1:n]
    topn_metrics = metrics[topn_positions]
    return topn_positions, topn_metrics
end

"""
    analogy_words(wv, pos, neg, n=5)

Return the top `n` words computed by analogy similarity between
positive words `pos` and negaive words `neg`. from the WordVectors `wv`.
"""
function analogy_words(wv::WordVectors, pos, neg, n=5)
    indx, metr = analogy(wv, pos, neg, n)
    return vocabulary(wv)[indx]
end

"""
    wordvectors(filename [,type=Float64][; kind=:text, header=false, normalize=true, vocabulary=nothing])

Generate a WordVectors type object from a file.

# Arguments
  * `filename::AbstractString` the embeddings file name
  * `type::Type` type of the embedding vector elements; default `Float64`

# Keyword arguments
  * `kind::Symbol` specifies whether the embeddings file is textual (`:text`)
or binary (`:binary`); default `:text`
  * `header::Union{Nothing, Bool}` in text embeddings files specifies
  whether the file contains a header i.e. number of lines, columns or not.
  If the header is `nothing`, the loader will attempt to autodetect the
  presence of a header; default `nothing`
  * `normalize:Bool` specifies whether to normalize the embedding vectors
i.e. return unit vectors; default true
  * `vocabulary::Union{Nothing, AbstractString}` path to the vocabulary
file generated with `vocab_count` (needed for binary embeddings);
default `nothing`
  * `load_bias::Bool` specifies whether to load the bias term or not
if using binary embedding files; default `false`
"""
function wordvectors(filename::AbstractString,
                     ::Type{T};
                     kind::Symbol=:text,
                     header::Union{Nothing, Bool}=nothing,
                     normalize::Bool=true,
                     vocabulary::Union{Nothing, AbstractString}=nothing,
                     load_bias::Bool=false
                    ) where T <: Real
    if kind == :binary
        vocabulary == nothing &&
            @error "A GloVe-generated vocabulary file path must be specified."
        return _from_binary(T, filename, vocabulary, load_bias, normalize)
    elseif kind == :text
        header == nothing && (header = autodetect_header(filename))
        if header
            return _from_text_header(T, filename, normalize)
        else
            return _from_text(T, filename, normalize)
        end
    else
        throw(ArgumentError("Unknown embedding file kind $(kind)"))
    end
end

"""
    autodetect_header(filename)

Function that attempts at autodetecting the presence of a header in a
GloVe embeddings file in a text format. If the function fails to detect
a header, a `false` value is returned.

Note: The function explicitly expects a text format for the embeddings.
The behaviour is undetermined for binary formats and no attempt at
detecting the file format is done.
"""
function autodetect_header(filename::AbstractString)
    open(filename) do f
        tokens1 = split(strip(readline(f)), ' ')
        tokens2 = split(strip(readline(f)), ' ')
        try
            # Attempt to parse the first line as header,
            # the second as an embedding
            @assert length(tokens1) == 2
            vocab_size = parse(Int, tokens1[1])
            vector_size = parse(Int, tokens1[2])
            vector = map(x-> parse(Float64, x), tokens2[2:end])
            @assert (vector_size == 2*(length(vector)+1) ||  # `model` is 0
                     vector_size == length(vector))          # `model` is 1 or 2
            return true
        catch
            return false
        end
    end
end

wordvectors(filename::AbstractString;
            kind::Symbol=:text,
            header::Bool=false,
            normalize::Bool=true,
            vocabulary::Union{Nothing, AbstractString}=nothing) =
    wordvectors(filename, Float64, kind=kind, header=header,
                normalize=normalize, vocabulary=vocabulary)

# Generate a WordVectors object from binary file
function _from_binary(::Type{T},
                      filename::AbstractString,
                      vocabulary::AbstractString,
                      load_bias::Bool=false,
                      normalize::Bool=true) where T<:Real
    vocab = Vector{String}(undef, 0)
    open(vocabulary) do f
        for line in eachline(f)
            word, count = split(line, ' ')
            push!(vocab, word)
        end
    end
    vocab_size = length(vocab)
    open(filename) do f
        nb = 0  # byte count
        while !eof(f)
            read(f,1)
            nb+= 1
        end
        vector_size = Int(nb / (2 * 8 * vocab_size))
        vector_mask = 1 : (vector_size - ifelse(load_bias, 0, 1))
        vectors = Matrix{T}(undef, length(vector_mask), vocab_size)
        binary_length = sizeof(Float64) * vector_size
        seekstart(f)
        for i in 1:vocab_size
            vector = reinterpret(Float64, read(f, binary_length))[vector_mask]
            if normalize
                vector = vector ./ norm(vector)  # unit vector
            end
            vectors[:,i] = T.(vector)
        end
        return WordVectors(vocab, vectors)
    end
end

# Generate a WordVectors object from text file
function _from_text_header(::Type{T},
                           filename::AbstractString,
                           normalize::Bool=true) where T<:Real
    open(filename) do f
        header = strip(readline(f))
        vocab_size, vector_size = map(x -> parse(Int, x), split(header, ' '))
        vocab = Vector{String}(undef, vocab_size+1)  # +1 for <unk> term
        vectors = Matrix{T}(undef, vector_size, vocab_size+1)
        i = 1
        for line in eachline(f)
            line = strip(line)
            parts = split(line, ' ')
            word = parts[1]
            vector = map(x-> parse(T, x), parts[2:end])
            if normalize
                vector = vector ./ norm(vector)  # unit vector
            end
            vocab[i] = word
            vectors[:, i] = vector
            i+= 1
        end
    return WordVectors(vocab, vectors)
    end
end

# Generate a WordVectors object from text file
function _from_text(::Type{T},
                    filename::AbstractString,
                    normalize::Bool=true) where T<:Real
    open(filename) do f
        vocab = Vector{String}()
        vectors = Vector{Vector{T}}()
        for (i, line) in enumerate(readlines(f))
            line = strip(line)
            parts = split(line, ' ')
            word = parts[1]
            vector = map(x-> parse(T, x), parts[2:end])
            if normalize
                vector = vector ./ norm(vector)  # unit vector
            end
            push!(vocab, word)
            push!(vectors, vector)
        end
        return WordVectors(vocab, reduce(hcat, vectors))
    end
end


text:
MOI.get(::Optimizer, ::MOI.SolverName) = "SeaPearl Solver"
MOI.get(model::Optimizer, ::MOI.TerminationStatus) = model.terminationStatus
MOI.get(model::Optimizer, ::MOI.PrimalStatus) = model.primalStatus

function MOI.get(optimizer::SeaPearl.Optimizer, ::MathOptInterface.VariablePrimal, vi::MathOptInterface.VariableIndex)
    id = vi.value
    last(optimizer.cpmodel.statistics.solutions)[string(id)]
end

"""
    MOI.is_empty(model::Optimizer)

Return a boolean saying if the model is empty or not.
"""
function MOI.is_empty(model::Optimizer)
    return isempty(model.cpmodel.variables) && isempty(model.cpmodel.constraints)
end

"""
    MOI.empty!(model::Optimizer)

Empty a given Optimizer.
"""
function MOI.empty!(model::Optimizer)
    # empty the cpmodel
    Base.empty!(model.cpmodel)
    model.variableselection = MinDomainVariableSelection()
    model.terminationStatus = MOI.OPTIMIZE_NOT_CALLED
    Base.empty!(model.moimodel.variables)
    Base.empty!(model.moimodel.constraints)
    Base.empty!(model.moimodel.affines)
    model.moimodel.objective_identifier = nothing

    # do not empty the options atm
end

MOI.supports(::Optimizer, ::MOI.RawParameter) = true
MOI.supports(::Optimizer, ::SeaPearl.MOIVariableSelectionAttribute) = true
MOI.supports(::Optimizer, ::SeaPearl.MOIValueSelectionAttribute) = true
MOI.supports(::Optimizer, ::MOI.VariablePrimal) = true

MOI.Utilities.map_indices(::Function, x::AbstractVariableSelection) = x
MOI.Utilities.map_indices(::Function, x::ValueSelection) = x

function MOI.set(model::Optimizer, p::MOI.RawParameter, value)
    model.options[p.name] = value
end

function MOI.set(model::Optimizer, ::SeaPearl.MOIVariableSelectionAttribute, heuristic::AbstractVariableSelection)
    model.variableselection = heuristic
end

function MOI.set(model::Optimizer, ::MOIValueSelectionAttribute, valueselection::ValueSelection)
    model.valueselection = valueselection
end

MOI.Utilities.supports_default_copy_to(model::Optimizer, copy_names::Bool) = !copy_names

function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike; kwargs...)
    return MOI.Utilities.automatic_copy_to(dest, src; kwargs...)
end


text:
@recipe function f(id::PerronFrobenius.InvariantDistribution)
    seriestype  :=  :bar
    xlabel --> "state # (i)"
    ylabel --> "invariant density"
    linecolor --> :black
    fillcolor --> :black
    fillalpha --> 0.5
    legend --> :none
    id.dist
end


text:
function Sinh()
    return Sinh(())
end

function value(obj::Sinh, arg0::Derivative)
    return jcall(obj, "value", Derivative, (Derivative,), arg0)
end

function value(obj::Sinh, arg0::DerivativeStructure)
    return jcall(obj, "value", DerivativeStructure, (DerivativeStructure,), arg0)
end

function value(obj::Sinh, arg0::jdouble)
    return jcall(obj, "value", jdouble, (jdouble,), arg0)
end


text:
module TestGroupBy
using AbstractTables
using Base.Test

include("mytable.jl")

name = ["Niamh", "Roger", "Genevieve", "Aiden"]
age = [27, 63, 26, 17]
eye_color = ["green", "brown", "brown", "blue"]

tbl = MyTable(
    name = NullableArray(name),
    age = NullableArray(age),
    eye_color = NullableArray(eye_color)
)
_tbl = copy(tbl)

group_indices = Dict{Any, Vector{Int}}([
    ("brown",false) => [3],
    ("green",true)  => [1],
    ("brown",true)  => [2],
    ("blue",false)  => [4],
])

qrya = @query groupby(tbl, eye_color, age > 26)
qryb = @query tbl |> groupby(eye_color, age > 26)

tbl2a = collect(qrya)
tbl2b = collect(qryb)

@test isequal(tbl, _tbl)
@test isequal(tbl2a, tbl2b)
for (group_level, indices) in zip(
    [
        (Nullable("brown"), Nullable(false)),
        (Nullable("green"), Nullable(true)),
        (Nullable("brown"), Nullable(true)),
        (Nullable("blue"), Nullable(false))
    ],
    [
        [3],
        [1],
        [2],
        [4]
    ]
)
    @test isequal(tbl2a.group_indices[group_level], indices)
end

@test isequal(tbl2a.groupbys, qrya.graph.args)

end


text:
# Taken from https://github.com/JuliaMath/KahanSummation.jl/pull/7/files
#   and make the follwoing changes: (1) Comment out the first two lines
#   ("__precom.." and "module...") and comment out all the lines after Line
#   100. (2) Update the function name `mapreduce_single` to `mapreduce_first`.
# Use `include("KahanSummation_patch.jl")` to load the patch.

# This file contains code that was formerly a part of Julia.
# License is MIT: https://julialang.org/license

# __precompile__(true)
# module KahanSummation

if VERSION >= v"0.7.0-DEV.3000" # TODO: More specific bound
    if isdefined(Base, :sum_kbn) # Deprecated
        import Base: sum_kbn, cumsum_kbn
    else
        export sum_kbn, cumsum_kbn
    end
end
if isdefined(Base, Symbol("@default_eltype"))
    using Base: @default_eltype
else
    macro default_eltype(itr)
        quote
            Core.Inference.return_type(first, Tuple{$(esc(itr))})
        end
    end
end
if isdefined(Base, :promote_sys_size_add)
    using Base: promote_sys_size_add
else
    promote_sys_size_add(x::T) where {T} = Base.r_promote(+, zero(T)::T)
end

"""
    TwicePrecisionN{T}
Represents an extended precision number as `x.hi - x.nlo`.
We store the lower order component as the negation to avoid problems when `x.hi == -0.0`.
"""
struct TwicePrecisionN{T}
    hi::T
    nlo::T
end

@inline function plus_kbn(x::T, y::T) where {T}
    hi = x + y
    nlo = abs(x) > abs(y) ? (hi - x ) - y : (hi - y) - x
    TwicePrecisionN(hi, nlo)
end
@inline function plus_kbn(x::T, y::TwicePrecisionN{T}) where {T}
    hi = x + y.hi
    if abs(x) > abs(y.hi)
        nlo = ((hi - x) - y.hi) + y.nlo
    else
        nlo = ((hi - y.hi) - x) + y.nlo
    end
    TwicePrecisionN(hi, nlo)
end
@inline plus_kbn(x::TwicePrecisionN{T}, y::T) where {T} = plus_kbn(y, x)

@inline function plus_kbn(x::TwicePrecisionN{T}, y::TwicePrecisionN{T}) where {T}
    hi = x.hi + y.hi
    if abs(x.hi) > abs(y.hi)
        nlo = (((hi - x.hi) - y.hi) + y.nlo) + x.nlo
    else
        nlo = (((hi - y.hi) - x.hi) + x.nlo) + y.nlo
    end
    TwicePrecisionN(hi, nlo)
end

Base.convert(::Type{TwicePrecisionN{T}}, x::Number) where {T} =
    TwicePrecisionN{T}(convert(T, x), zero(T))
Base.convert(::Type{T}, x::TwicePrecisionN) where {T} =
    convert(T, x.hi - x.nlo)

@static if VERSION >= v"0.7.0-"
    Base.mapreduce_empty(f, ::typeof(plus_kbn), T) = TwicePrecisionN(zero(T),zero(T))
    Base.mapreduce_empty(::typeof(identity), ::typeof(plus_kbn), T) = TwicePrecisionN(zero(T),zero(T)) # disambiguate
    Base.mapreduce_first(f, ::typeof(plus_kbn), x) = TwicePrecisionN(x, zero(x))
else
    Base.r_promote_type(::typeof(plus_kbn), ::Type{T}) where {T} =
        TwicePrecisionN{T}
    Base.mr_empty(f, ::typeof(plus_kbn), T) = TwicePrecisionN(zero(T),zero(T))
end

singleprec(x::TwicePrecisionN{T}) where {T} = convert(T, x)

"""
    sum_kbn([f,] A)
Return the sum of all elements of `A`, using the Kahan-Babuska-Neumaier compensated
summation algorithm for additional accuracy.
"""
sum_kbn(f, X) = singleprec(mapreduce(f, plus_kbn, X))
sum_kbn(X) = sum_kbn(identity, X)

#=

"""
    cumsum_kbn(A, dim::Integer)
@@ -85,32 +157,4 @@ function cumsum_kbn(v::AbstractVector{T}) where T<:AbstractFloat
    return r
end

"""
    sum_kbn(A)
Return the sum of all elements of `A`, using the Kahan-Babuska-Neumaier compensated
summation algorithm for additional accuracy.
"""
function sum_kbn(A)
    T = @default_eltype(typeof(A))
    c = promote_sys_size_add(zero(T)::T)
    i = start(A)
    if done(A, i)
        return c
    end
    Ai, i = next(A, i)
    s = Ai - c
    while !(done(A, i))
        Ai, i = next(A, i)
        t = s + Ai
        if abs(s) >= abs(Ai)
            c -= ((s-t) + Ai)
        else
            c -= ((Ai-t) + s)
        end
        s = t
    end
    s - c
end

end # module
=#


text:
#' ## Functions
#'
#'
#' In the last notebook, we talked about modeling data with functions. A **function** is one of the most fundamental concepts in computing (and also in mathematics).
#'
#' A function is a piece of a program that receives **input arguments**, processes them by doing certain calculations on them, and returns **outputs**.
#'
#' For example, we might have a function `g` that takes a number as an input and returns the square of that number as an output. How can we define this function `g` on a computer? Julia gives us a few different ways to do this.
#' ### Defining functions
#' Firstly, we could write `g` as follows:

#+ results="hidden"
g(x) = x^2
#+
#'
g(2), g(3.5) # This is a tuple.
#'
g(3.0)
#' Alternatively, we could declare this function using the `function` and `end` keywords:
#+ results="hidden"
function g1(x)
    x^2
end
#+
g1(2), g1(3.5)
#' The third way we could have declared this function is as an "anonymous" or
#' "lambda" function. "Anonymous" functions are functions that truly don't need
#' names! For example, we could have declared a function that squares its input
#' as
(x->x^2)(2)
#' Now that we've done that, we can't access the function `x -> x^2` again
#' because we have no name to call! That seems a little silly, doesn't it?
#'
#' Actually, there are times where functions without names are especially
#' handy.  Most commonly they're used arguments to "higher-order" functions:
#' these are functions which take _other functions_ as arguments. For example:
map(sqrt, [1, 2, 3])
#'
map(x->x^2, [1, 2, 3])
#' Of course, you can also assign anonymous function a name, but in general
#' it's better style to use the non-anonymous equivalents above if you actually
#' want to use a name.
g2 = x -> x^2
g2(3.5), g2("I ♡ Julia ") # Use \heartsuit + TAB to get the ♡ character

#' ## An important sigmoidal function
#' A particular function that is used a lot in machine learning is a so-called
#' "sigmoidal" function (meaning a function that is S-shaped, i.e. the graph
#' of the function looks like an `S`).
#'
#' $$\sigma(x) = \frac{1}{1 + e^{-x}}$$
#'
#' The sigmoid function that we will use is given the name $\sigma$, and is
#' defined by the following mathematical expression:
#'
#' $$\sigma(x) := \frac{1}{1 + \exp(-x)}$$

#' #### Exercise 1
#'
#' Use the first syntax given above to define the function `σ` in Julia.
#' Note that Julia actually allows us to use the symbol σ as a variable name!
#' To do so, type `\sigma<TAB>` in the code cell.

#' #### Solution
# Short form
σ(x) = 1/(1 + exp(-x))
σ(1)
#+ results="hidden"
# Long form
function sigma(x)
    1/(1 + exp(-x))
end
#+
sigma(1)

#' ## Plotting functions

#' Let's draw the function σ to see what it looks like. Throughout this course,
#' we'll use the Julia package `Plots.jl` for all of the graphics. This package
#' provides a flexible syntax for plotting, in which options to change
#' attributes like the width of the lines used in the figure are given as named
#' keyword arguments.
#'
#' In addition, it allows us to use different "backends", which are the other
#' libraries that actually carry out the plotting following the instructions
#' from `Plots.jl`.
using Plots
#gr()
# plotly interferes mathjax, tpl file needs to be fixed
# see templates/julia_html.tpl
plotly()

pl = plot(sigma, -5, 5, legend=(0.1,0.9), label="σ")
# add horizontal lines at 0 and 1, with dashed style and linewidth 3
hline!(pl, [0,1], ls=:dash, lw=3, label=false)
# add a vertical line at 0
vline!(pl, [0], ls=:dash, lw=3, label=false)

#' We can think of $\sigma$ as a smooth version of a step or threshold function
#' (often called a "Heaviside" function). To see this, let's modify the
#' steepness of the jump in $\sigma$ and compare it to the Heaviside function;
#' we'll see how all this works in more detail later:
#+ results="hidden"
function heaviside(x)
    x < 0 ? 0.0 : 1.0
end
#+

#' Try manipulating the value of `w` between 0 to :
w = 10.0
plot(x->sigma(w*x), -5, 5, label="σ", lw=2, legend=(0.1,0.9))
plot!(heaviside, ls=:dash, label="step")
# This particular function takes any real number as input, and gives an output
# between $0$ and $1$. It is continuous and smooth.
#' Try a smaller $w$.
w = 1.0
plot(x->sigma(w*x), -5, 5, label="σ", lw=2, legend=(0.1,0.9))
plot!(heaviside, ls=:dash, label="step")
#' #### Exercise 2
#' Declare the sigmoid function above as an anonymous function with a different
#' name.
σx = x -> 1/(1+exp(-x))
σx(1)
#' ### Mutating functions: `...!`
#' To generate our plot of σ above, we used some functions that end with `!`.
#'  What does a `!` at the end of a function name mean in Julia?
#'
#' Functions that change or modify their inputs are called
#' **mutating functions**. But wait, don't all functions do that?
#' No, actually. Functions typically take *inputs* and use those *inputs* to
#' generate *outputs*, but the inputs themselves usually don't actually get
#' changed by a function. For example, copy and execute the following code:
v1 = [9, 4, 7, 11]
v2 = sort(v1)
#' `v2` is a sorted version of `v1`, but after calling `sort`, `v1` is still
#' unsorted.
v1
#+ results="hidden"
sort!(v1)
#' Look at the values in `v1` now!
v1
#' This time, the original vector itself was changed (mutated), and is now
#' sorted. Unlike `sort`, `sort!` is a mutating function. Did the `!` make
#' `sort!` mutating? Well, no, not really. In Julia, `!` indicates mutating
#' functions by convention. When the author of `sort!` wrote `sort!`, they added
#' a `!` to let you to know that `sort!` is mutating, but the `!` isn't what
#' makes a function mutating or non-mutating in the first place.

#' #### Exercise
#' Some of our plotting commands end with `!`. Copy and execute the following code:
r = -5:0.1:5
g(x) = x^2
h(x) = x^3
plot(r, g, label="g")
plot!(r, h, label="h")
#' Then change the code slightly to remove the `!` after `plot!(r, h)`.
#' How does this change your output? What do you think it means to add `!`
#' after plotting commands?
#' #### Solution
r = -5:0.1:5
g(x) = x^2
h(x) = x^3
p1 = plot(r, g, label="g")
#'
p2 = plot(r, h, label="h")
#' ```julia
#' plot(r, g)
#' plot!(r, h)
#' ```
#' creates an overlay of `g` and `h`, whereas
#'
#' ```julia
#' plot(r, g)
#' plot(r, h)
#' ```
#'
#' creates one plot for `g` and one for `h`.
#'
#' When we add a `!` after a plotting command, we are mutating or updating an
#' *existing* plot.

#' ## Pointwise application of functions, `f.(x, y)` and `x .+ y` - broadcasting
#' We saw in a previous notebook that we needed to add `.` after the names of
#' some functions, as in
#'
#' ```julia
#' green_amount = mean(Float64.(green.(apple)))
#' ```
#'
#' What are those extra `.`s really doing?
#'
#' When we add a `.` after a function's name, we are telling Julia that we want
#' to "**broadcast**" that function over the inputs passed to the function.
#' This means that we want to apply that function *element-wise* over the
#' inputs; in other words, it will apply the function to each element of the
#' input, and return an array with the newly-calculated values.

#' For example, copy and execute the following code:
g.(r)
#' Since the function `g` squares it's input, this squares all the elements of
#' the range `r`.
#' What happens if instead we just call `g` on `r` via
# g(r)

#' You should see an error message after calling `g(r)`, which says that Julia
#' cannot multiply two vectors. When we call `g(r)`, we ask Julia to multiply
#' `r` by `r`. When we call `g.(r)`, we ask Julia to multiply *each element*
#' in `r` by itself.

#' Try this and see what happens.
x = [1 2 3; 4 5 6]
f(t) = sqrt(t) + 5
f.(x)

#' #### Exercise 3
#'
#' Copy and execute the following code to get the type of the object
#' `numbers = [1, 2, "three", 4.0]`:
#+ results="hidden"
numbers = [1, 2, "three", 4.0]
#' What is the type of `numbers`?
typeof(numbers)
#' Its type is `Array{Any,1}`. This means that `numbers` is a
#' 1-dimensional array that has elements of abstract type `Any`. We can think
#' of the `Any` type as a classification that includes *all* concrete types,
#' for example, `Float64`, `Int32`, `Bool`, `String`, `Char`, etc.

#' #### Exercise 4
#'
#' Broadcast `typeof` over `numbers` to see what the types of the elements
#' stored inside `numbers` are.
typeof.(numbers)

#' Note: Alternatively, we could have looked at the `typeof` each element in
#' the array `numbers` using a `for` loop:
for n in numbers
    println(typeof(n))
end

#' #### Exercise 5
#'
#' Write a `for` loop that applies `g` to each of the elements of `r` and
#' prints the results. Verify that the numbers printed by this `for` loop are
#' equal to the entries of the output of `g.(r)`.
r = -5:0.1:5
for x in r
    #println(g(x)) # uncomment to print
end

#' We can check this via
broadcasted_g = g.(r)
for (x, bg) = zip(r, broadcasted_g)
    # println(g(x) == bg) # uncomment to print
end

#' A better way
answers_match = true
for (x, bg) = zip(r, broadcasted_g)
    if (g(x) != bg)
        answers_match = false
    end
end
answers_match ? println("Our `for` loop worked!") : println("Our `for` loop does not work")

#' #### Exercise 6
#'
#' Define a range `xs` between -5 and 5 with steps of 0.5.
#' Apply the $\sigma$ function pointwise to this range and define `ys` as the
#' result.
#' What does the result look like? Plot these as points and join them with
#' lines.
#'
#' Make the plot interactive where you can vary the step size. Fix the range of
#' the plot in the `x` and `y` directions using the functions `xlims!` and
#' `ylims!`.
#'
#' Learn about plot optinos by `plotattr(:Series)` and `plotattr("markershape")`
xs = -5:0.5:5
ys = sigma.(xs)
pl = plot(xs, ys; legend=(0.1,0.9), markershape=:circle)

#' Try manipulating `stepsize` between 0 and 1:
stepsize = 0.5
xs = -5.0:stepsize:5.0
ys = σ.(xs)

scatter(xs, ys)
plot!(xs, ys)

xlims!(-5, 5)
ylims!(0, 1)


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to the Euler angle and axis.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

################################################################################
#                                  Operations
################################################################################

"""
    *(av₂::EulerAngleAxis{T1}, av₁::EulerAngleAxis{T2}) where {T1,T2}

Compute the composed rotation of `av₁ -> av₂`.

The rotation will be represented by a Euler angle and axis (see
[`EulerAngleAxis`](@ref)). By convention, the output angle will always be in the
range `[0, π] rad`.

Notice that the vector representing the axis in `av₁` and `av₂` must be unitary.
This function neither verifies this nor normalizes the vector.

# Examples

```jldoctest
julia> av1 = EulerAngleAxis(deg2rad(45), [sqrt(2)/2, sqrt(2)/2, 0])
EulerAngleAxis{Float64}:
  Euler angle : 0.785398 rad  (45.0°)
  Euler axis  : [0.707107, 0.707107, 0.0]

julia> av2 = EulerAngleAxis(deg2rad(22.5), [sqrt(2)/2, sqrt(2)/2, 0])
EulerAngleAxis{Float64}:
  Euler angle : 0.392699 rad  (22.5°)
  Euler axis  : [0.707107, 0.707107, 0.0]

julia> av1 * av2
EulerAngleAxis{Float64}:
  Euler angle : 1.1781 rad  (67.5°)
  Euler axis  : [0.707107, 0.707107, 0.0]
```
"""
function *(av₂::EulerAngleAxis{T1}, av₁::EulerAngleAxis{T2}) where {T1, T2}
    # Auxiliary variables.
    sθ₁o2, cθ₁o2 = sincos(av₁.a / 2)
    sθ₂o2, cθ₂o2 = sincos(av₂.a / 2)

    v₁ = av₁.v
    v₂ = av₂.v

    # Compute `cos(θ/2)` in which `θ` is the new Euler angle.
    cθo2 = cθ₁o2 * cθ₂o2 - sθ₁o2 * sθ₂o2 * dot(v₁, v₂)

    T = promote_type(T1, T2, typeof(sθ₁o2), typeof(sθ₂o2), typeof(cθo2))

    if abs(cθo2) >= 1 - eps()
        # In this case, the rotation is the identity.
        return EulerAngleAxis(T(0), SVector{3,T}(0, 0, 0))
    else
        # Compute `sin(θ/2)` in which `θ` is the new Euler angle.
        sθo2 = √(1 - cθo2 * cθo2)

        # Compute the θ angle between [0, 2π].
        θ = 2acos(cθo2)

        # Keep the angle between [0, π].
        s = +1
        if θ > π
            θ = T(2π) - θ
            s = -1
        end

        v = s * (sθ₁o2 * cθ₂o2 * v₁ + cθ₁o2 * sθ₂o2 * v₂ + sθ₁o2 * sθ₂o2 * (v₁ × v₂)) / sθo2

        return EulerAngleAxis(θ, v)
    end
end

"""
    inv(av::EulerAngleAxis)

Compute the inverse rotation of the Euler angle and axis `av`.

The Euler angle returned by this function will always be in the interval `[0, π]
rad`.

# Examples

```jldoctest
julia> av = EulerAngleAxis(deg2rad(20), [sqrt(2) / 2, 0, sqrt(2) / 2])
EulerAngleAxis{Float64}:
  Euler angle : 0.349066 rad  (20.0°)
  Euler axis  : [0.707107, 0.0, 0.707107]

julia> inv(av)
EulerAngleAxis{Float64}:
  Euler angle : 0.349066 rad  (20.0°)
  Euler axis  : [-0.707107, -0.0, -0.707107]

julia> av = EulerAngleAxis(deg2rad(-20), [sqrt(2) / 2, 0, sqrt(2) / 2])
EulerAngleAxis{Float64}:
  Euler angle : -0.349066 rad  (-20.0°)
  Euler axis  : [0.707107, 0.0, 0.707107]

julia> inv(av)
EulerAngleAxis{Float64}:
  Euler angle : 0.349066 rad  (20.0°)
  Euler axis  : [0.707107, 0.0, 0.707107]
```
"""
@inline function inv(av::EulerAngleAxis{T}) where T<:Number
    # Make sure that the Euler angle is always in the inverval [0,π]
    s = -1
    θ = mod(av.a, T(2π))

    if θ > π
        s = 1
        θ = T(2π) - θ
    end

    return EulerAngleAxis(θ, s*av.v)
end

################################################################################
#                                      IO
################################################################################

function show(io::IO, av::EulerAngleAxis{T}) where T
    # Get if `io` request a compact printing, defaulting to true.
    compact_printing = get(io, :compact, true)

    # Convert the values using `print` and compact printing.
    θ_str  = sprint(print, av.a;    context = :compact => compact_printing)
    v₁_str = sprint(print, av.v[1]; context = :compact => compact_printing)
    v₂_str = sprint(print, av.v[2]; context = :compact => compact_printing)
    v₃_str = sprint(print, av.v[3]; context = :compact => compact_printing)

    print(io, "EulerAngleAxis{$T}: θ = " * θ_str * " rad, v = [" *
        v₁_str * ", " * v₂_str * ", " * v₃_str * "]")

    return nothing
end

function show(io::IO, mime::MIME"text/plain", av::EulerAngleAxis{T}) where T
    # Check if the user wants colors.
    color = get(io, :color, false)::Bool

    # Check if the user wants compact printing, defaulting to `true`.
    compact_printing = get(io, :compact, true)::Bool

    # Assemble the context.
    context = IOContext(io, :color => color, :compact => compact_printing)

    g = color ? string(_g) : ""
    y = color ? string(_y) : ""
    d = color ? string(_d) : ""

    θ_str  = sprint(print, av.a; context = context)
    θd_str = sprint(print, rad2deg(av.a); context = context)
    v₁_str = sprint(print, av.v[1]; context = context)
    v₂_str = sprint(print, av.v[2]; context = context)
    v₃_str = sprint(print, av.v[3]; context = context)

    println(io, "EulerAngleAxis{$T}:")
    print(  io, g)
    print(  io, "  Euler angle : ")
    print(  io, d)
    print(  io, θ_str * " rad  ")
    println(io, "(" * θd_str * "°)")
    print(  io, y)
    print(  io, "  Euler axis  : ")
    print(  io, d)
    print(io, "[" * v₁_str * ", " * v₂_str * ", " * v₃_str * "]")

    return nothing
end


text:
addSaving_! =
Dict{String,Function}(
    "RAM"=>addSavingRAM_!,
    "CSV"=>addSavingCSV_!,
)

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia