text:
function testNHDensityCurrent()
close all
clear all

% Physical parameters
Param=PhysParameters();

nx=60;
ny=2;
lx=2*25600;
ly=2000;
x0=-25600;
y0=0;
Boundary.WE='Period';
Boundary.BT='Period';
Param.hS='';
Param.Grid=CartGrid(nx,ny,lx,ly,x0,y0,OrientFaceCart,Boundary,Param);

Param.H=6400;
nz=40;
zP=zeros(nz,1);
z=zeros(nz+1,1);
dz=Param.H/nz;
zP(1)=dz/2;
for i=2:nz
  zP(i)=zP(i-1)+dz;
end
for i=2:nz+1
  z(i)=z(i-1)+dz;
end

Param.Grid.nz=nz;
Param.Grid.z=z;
Param.Grid.zP=zP;
Param.Grid.dz=dz;

% Model
Param.ModelType='Curl';
%Param.ModelType='Div';
Param.Thermo='';
Param.Source=false;
Param.Damping=false;
Param.Flat=false;
Param.Coriolis=false;
Param.Buoyancy=true;
Param.xC0=0;
Param.zC0=3000;
Param.xrC0=4000;
Param.zrC0=2000;
Param.T0=300;
Param.DeltaT=-15.0;
Param.uMax=0;
Param.vMax=0;
Param.NBr=1.e-2;
Param.Equation='Compressible';
Param.TopoS='';
Param.EndTime=900;

% Discretization
OrdPoly=4;
OrdPolyZ=1;
[CG,Param]=Discretization(OrdPoly,OrdPolyZ,@JacobiDG3,Param);
LRef=11*1.e5;
Param.HyperVisc=true;
Param.HyperDCurl=2.e5; %1.e14*(dx/LRef)^3.2;
Param.HyperDGrad=2.e5;
Param.HyperDDiv=2.2e5; %2.e2;

% Initial conditions
Param.NumV=5;
U=zeros(CG.NumG,nz,Param.NumV);
Param.ProfRho='DensityCurrent';
Param.ProfTheta='DensityCurrent';
Param.ProfVel='Const';
Param.RhoPos=1;
Param.uPos=2;
Param.vPos=3;
Param.wPos=4;
Param.ThPos=5;
U(:,:,Param.RhoPos)=Project(@fRho,CG,Param);
[U(:,:,Param.uPos),U(:,:,Param.vPos)]=ProjectVec(@fVel,CG,Param);
U(:,:,Param.ThPos)=Project(@fTheta,CG,Param).*U(:,:,Param.RhoPos);

% Output
Slice.Type='XZ';
Slice.y=y0+ly/ny/2;
Slice.iy=1;
Param.SliceXZ=Slice;
Param.SliceXY.Type='XY';
Param.SliceXY.iz=1;
Param.fig=1;
Param.vtk=0;
Param.vtkFileName='DensityCurrent';
vtkGrid=vtkCGGrid(CG,@TransCart,@Topo,Param);
Param.cNames(1).s='Rho';
Param.cNames(2).s='u';
Param.cNames(3).s='w';
Param.cNames(4).s='Th';
cOut=zeros(CG.NumG,nz,4);

%IntMethod='Rosenbrock';
IntMethod='RungeKutta';
dtau=.2;
Param.EndTime=900;
PrintTime=30;
%v = VideoWriter ('Galewsky.avi');
%open (v);
Param.RK=RungeKuttaMethod('RK4');
CFL=0.125;
time=0;
nIter=Param.EndTime/dtau;
PrintInt=PrintTime/dtau;
cOut(:,:,1)=U(:,:,Param.RhoPos);
cOut(:,:,2)=U(:,:,Param.uPos);
cOut(:,1,3)=0.5*U(:,1,Param.wPos);
cOut(:,2:nz-1,3)=0.5*(U(:,1:nz-2,Param.wPos)+U(:,2:nz-1,Param.wPos));
cOut(:,nz,3)=0.5*U(:,nz-1,Param.wPos);
cOut(:,:,4)=U(:,:,Param.ThPos)./U(:,:,Param.RhoPos);
Param.vtk=vtkCG(cOut,CG,Param,vtkGrid,Param.vtk);
switch IntMethod
  case 'RungeKutta'
    for i=1:nIter
      i
      U=RungeKuttaExplicit(U,dtau,@FcnNHCurlVec,CG,Param);
      time=time+dtau;
      if mod(i,PrintInt)==0
        cOut(:,:,1)=U(:,:,Param.RhoPos);
        cOut(:,:,2)=U(:,:,Param.uPos);
        cOut(:,1,3)=0.5*U(:,1,Param.wPos);
        cOut(:,2:nz-1,3)=0.5*(U(:,1:nz-2,Param.wPos)+U(:,2:nz-1,Param.wPos));
        cOut(:,nz,3)=0.5*U(:,nz-1,Param.wPos);
        cOut(:,:,4)=U(:,:,Param.ThPos)./U(:,:,Param.RhoPos);
        Param.vtk=vtkCG(cOut,CG,Param,vtkGrid,Param.vtk);
      end
    end
end
%fig=PlotCG(U(:,1),CG,@JacobiSphere2,Param,fig);
%close(v)
end


text:
using Pkg; Pkg.activate("."); using Revise;
using Plots; using IntensityScans; using IntensityScans.ThorlabsLTStage; using IntensityScans.TcpInstruments; using IntensityScans.Unitful;

params = ScanParameters(
    medium = Medium(),
    excitation = Excitation(),
    f0 = 15e6,
    hydrophone_id = :Onda_HGL0200_2322,
    preamp_id = :Onda_AH2020_1238_20dB,
)

info_xy = load("native_vsx_subaperature_focusing.jld2")
info_xy = compute_metrics(info_xy, params)

ext: jl


text:
using DelimitedFiles
@testset "Testing main functions" begin

options = ConformationSetup(0.00001,classicBP,true)
data = preprocessing("ins_15v2.pdb")

sol = conformation(data,options)

@test sol.number == 4

mol1 = MoleculeType([AtomType(0.0, 0.0, 0.0),
		     AtomType(-1.526, 0.0, 0.0),
		     AtomType(-2.034999999999995,1.4379999999999997,0.0),
		     AtomType(-1.6749999999999972,2.1060000000000003,-1.3239999999999965),
		     AtomType(-2.4299999999999944,3.4260000000000015,-1.4499999999999982),
		     AtomType(-1.9960000000000164,4.140999999999982,-2.726000000000016),
		     AtomType(-0.47700000000001186,4.283999999999978,-2.7380000000000426),
		     AtomType(-0.05900000000000982,5.365999999999996,-1.746000000000059),
		     AtomType(-0.6160000000000029,6.71299999999999,-2.198000000000081),
		     AtomType(-2.138000000000004,6.632999999999996,-2.283000000000105),
		     AtomType(-2.7060000000000017,6.282000000000035,-0.9110000000000911),
		     AtomType(-4.211000000000007,6.055000000000045,-1.0260000000000806),
		     AtomType(-4.863999999999999,7.281000000000017,-1.6600000000001363),
		     AtomType(-5.922122878027752,7.841170383964165,-0.7130897001703507),
		     AtomType(-6.0663693289671485,9.342642843087825,-0.9430775604685648)], 7.862256122994717e-28 )

@test sol.molecules[1]≈mol1

mol2 =  MoleculeType([AtomType( 0.0 , 0.0 ,  0.0 ),
		      AtomType( -1.526 , 0.0 ,  0.0 ),
		      AtomType( -2.034999999999995 , 1.4379999999999997 ,  0.0 ),
		      AtomType( -1.6749999999999972 , 2.1060000000000003 ,  -1.3239999999999965 ) ,
		      AtomType( -2.4299999999999944 , 3.4260000000000015 ,  -1.4499999999999982 ) ,
		      AtomType( -1.9960000000000164 , 4.140999999999982 ,  -2.726000000000016 ) ,
		      AtomType( -0.47700000000001186 , 4.283999999999978 ,  -2.7380000000000426 ),
		      AtomType( -0.05900000000000982 , 5.365999999999996 ,  -1.746000000000059 ) ,
		      AtomType( -0.6160000000000029 , 6.71299999999999 ,  -2.198000000000081 ) ,
		      AtomType( -2.138000000000004 , 6.632999999999996 ,  -2.283000000000105 ) ,
		      AtomType( -2.7060000000000017 , 6.282000000000035 ,  -0.9110000000000911 ),
		      AtomType( -4.211000000000007 , 6.055000000000045 ,  -1.0260000000000806 ) ,
		      AtomType( -4.863999999999999 , 7.281000000000017 ,  -1.6600000000001363 ) ,
		      AtomType( -6.302 , 6.948000000000005 ,  -2.0490000000001407 ) ,
		      AtomType( -7.122000000000005 , 6.687000000000051 ,  -0.7890000000001312 )],7.853878489881605e-28)

@test sol.molecules[2]≈mol2

mol3 = MoleculeType([AtomType( 0.0 , 0.0 ,  0.0 ) ,
		     AtomType( -1.526 , 0.0 ,  0.0 ) ,
		     AtomType( -2.034999999999995 , 1.4379999999999997 ,  0.0 ) ,
		     AtomType( -1.6749999999999972 , 2.1060000000000003 ,  1.3239999999999965 ) ,
		     AtomType( -2.4299999999999944 , 3.4260000000000015 ,  1.4499999999999982 ) ,
		     AtomType( -1.9960000000000164 , 4.140999999999982 ,  2.726000000000016 ) ,
		     AtomType( -0.47700000000001186 , 4.283999999999978 ,  2.7380000000000426 ),
		     AtomType( -0.05900000000000982 , 5.365999999999996 ,  1.746000000000059 ) ,
		     AtomType( -0.6160000000000029 , 6.71299999999999 ,  2.198000000000081 ) ,
		     AtomType( -2.138000000000004 , 6.632999999999996 ,  2.283000000000105 ) ,
		     AtomType( -2.7060000000000017 , 6.282000000000035 ,  0.9110000000000911 ),
		     AtomType( -4.211000000000007 , 6.055000000000045 ,  1.0260000000000806 ) ,
		     AtomType( -4.863999999999999 , 7.281000000000017 ,  1.6600000000001363 ) ,
		     AtomType( -6.302 , 6.948000000000005 ,  2.0490000000001407 ) ],7.852246476949892e-28)

@test sol.molecules[3]≈mol3

mol4 = MoleculeType([AtomType( 0.0 , 0.0 ,  0.0 ),
		     AtomType( -1.526 , 0.0 ,  0.0 ),
		     AtomType( -2.034999999999995 , 1.4379999999999997 ,  0.0 ),
		     AtomType( -1.6749999999999972 , 2.1060000000000003 ,  1.3239999999999965 ),
		     AtomType( -2.4299999999999944 , 3.4260000000000015 ,  1.4499999999999982 ),
		     AtomType( -1.9960000000000164 , 4.140999999999982 ,  2.726000000000016 ),
		     AtomType( -0.47700000000001186 , 4.283999999999978 ,  2.7380000000000426 ),
		     AtomType( -0.05900000000000982 , 5.365999999999996 ,  1.746000000000059 ),
		     AtomType( -0.6160000000000029 , 6.71299999999999 ,  2.198000000000081 ),
		     AtomType( -2.138000000000004 , 6.632999999999996 ,  2.283000000000105 ),
		     AtomType( -2.7060000000000017 , 6.282000000000035 ,  0.9110000000000911 ),
		     AtomType( -4.211000000000007 , 6.055000000000045 ,  1.0260000000000806 ),
		     AtomType( -4.863999999999999 , 7.281000000000017 ,  1.6600000000001363 ),
		     AtomType( -5.922122878027752 , 7.841170383964165 ,  0.7130897001703507 ),
		     AtomType( -7.059383306135178 , 6.834224913707539 ,  0.5687964619691643 )],7.863165562712991e-28)

@test sol.molecules[4]≈mol4

end


text:
module CacheVariables

using BSON
import Logging: @info

export @cache, cache

function _cachevars(ex::Expr)
    (ex.head === :(=)) && return Symbol[ex.args[1]]
    (ex.head === :block) && return collect(Iterators.flatten([
        _cachevars(exi) for exi in ex.args if isa(exi, Expr)
    ]))
    return Vector{Symbol}(undef, 0)
end

"""
    @cache path code [overwrite]

Cache results from running `code` using BSON file at `path`.
Load if the file exists; run and save if it does not.
Run and save either way if `overwrite` is true (default is false).

Tip: Use `begin...end` for `code` to cache blocks of code.

Caveat: The variable name `ans` is used for storing the final output,
so it is best to avoid using it as a variable name in `code`.

# Examples
```julia-repl
julia> @cache "test.bson" begin
         a = "a very time-consuming quantity to compute"
         b = "a very long simulation to run"
         100
       end
┌ Info: Saving to test.bson
│ a
└ b
100

julia> @cache "test.bson" begin
         a = "a very time-consuming quantity to compute"
         b = "a very long simulation to run"
         100
       end
┌ Info: Loading from test.bson
│ a
└ b
100

julia> @cache "test.bson" begin
         a = "a very time-consuming quantity to compute"
         b = "a very long simulation to run"
         100
       end true
┌ Info: Overwriting test.bson
│ a
└ b
100
```
"""
macro cache(path, ex::Expr, overwrite = false)
    vars = _cachevars(ex)
    vardesc = join(string.(vars), "\n")
    varkws = [:($(var) = $(var)) for var in vars]
    varlist = :($(varkws...),)
    vartuple = :($(vars...),)

    return quote
        if !ispath($(esc(path))) || $(esc(overwrite))
            _ans = $(esc(ex))
            _msg = ispath($(esc(path))) ? "Overwriting " : "Saving to "
            @info(string(_msg, $(esc(path)), "\n", $(vardesc)))
            mkpath(splitdir($(esc(path)))[1])
            bson($(esc(path)); $(esc(varlist))..., ans = _ans)
            _ans
        else
            @info(string("Loading from ", $(esc(path)), "\n", $(vardesc)))
            data = BSON.load($(esc(path)), @__MODULE__)
            $(esc(vartuple)) = getindex.(Ref(data), $vars)
            data[:ans]
        end
    end
end

"""
    cache(f, path; mod = @__MODULE__)

Cache output from running `f()` using BSON file at `path`.
Load if the file exists; run and save if it does not.
Use `mod` keyword argument to specify module.

Tip: Use `do...end` to cache output from a block of code.

# Examples
```julia-repl
julia> cache("test.bson") do
         a = "a very time-consuming quantity to compute"
         b = "a very long simulation to run"
         return (; a = a, b = b)
       end
[ Info: Saving to test.bson
(a = "a very time-consuming quantity to compute", b = "a very long simulation to run")

julia> cache("test.bson") do
         a = "a very time-consuming quantity to compute"
         b = "a very long simulation to run"
         return (; a = a, b = b)
       end
[ Info: Loading from test.bson
(a = "a very time-consuming quantity to compute", b = "a very long simulation to run")
```
"""
function cache(@nospecialize(f), path; mod = @__MODULE__)
    if !ispath(path)
        ans = f()
        @info(string("Saving to ", path, "\n"))
        mkpath(splitdir(path)[1])
        bson(path; ans = ans)
        return ans
    else
        @info(string("Loading from ", path, "\n"))
        data = BSON.load(path, mod)
        return data[:ans]
    end
end

end # module


text:
module PythonCallChainRules

using Requires

import FillArrays
import Adapt

struct PyAdaptor{T} end
Adapt.adapt_storage(to::PyAdaptor{T}, x::AbstractArray) where {T} = convert(Array, x)
Adapt.adapt_storage(to::PyAdaptor{T}, x::StridedArray) where {T} = x
Adapt.adapt_storage(to::PyAdaptor{T}, x::FillArrays.AbstractFill) where {T} = collect(x)

include("pytorch.jl")
include("jax.jl")

function __init__()
    @require CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" begin
        include("cuda.jl")
    end
end

end


text:
using Statistics: mean!, sum!

using Oceananigans.Utils: tupleit
using Oceananigans.Grids: regular_dimensions
using Oceananigans.Fields: condition_operand
import Oceananigans.Fields: Reduction

#####
##### Metric inference
#####

reduction_grid_metric(dims::Number) = reduction_grid_metric(tuple(dims))

reduction_grid_metric(dims) = dims === tuple(1)  ? Δx :
                              dims === tuple(2)  ? Δy :
                              dims === tuple(3)  ? Δz :
                              dims === (1, 2)    ? Az :
                              dims === (1, 3)    ? Ay :
                              dims === (2, 3)    ? Ax :
                              dims === (1, 2, 3) ? volume :
                              throw(ArgumentError("Cannot determine grid metric for reducing over dims = $dims"))

#####
##### Metric reductions
#####

struct Average end

function Reduction(avg::Average, field::AbstractField; condition = nothing, mask = 0, dims)
    dims = dims isa Colon ? (1, 2, 3) : tupleit(dims)
    dx = reduction_grid_metric(dims)

    if all(d in regular_dimensions(field.grid) for d in dims)
        # Dimensions being reduced are regular; just use mean!
        return Reduction(mean!, condition_operand(field, condition, mask); dims)
    else
        # Compute "size" (length, area, or volume) of averaging region
        metric = GridMetricOperation(location(field), dx, field.grid)
        L = sum(metric; condition, mask, dims)

        # Construct summand of the Average
        L⁻¹_field_dx = field * dx / L

        return Reduction(sum!, condition_operand(L⁻¹_field_dx, condition, mask), dims)
    end
end

"""
    Average(field::AbstractField; condition = nothing, mask = 0, dims=:)

Return `Reduction` representing a spatial average of `field` over `dims`.

Over regularly-spaced dimensions this is equivalent to a numerical `mean!`.

Over dimensions of variable spacing, `field` is multiplied by the
appropriate grid length, area or volume, and divided by the total
spatial extent of the interval.
"""
Average(field::AbstractField; condition = nothing, mask = 0, dims=:) = Reduction(Average(), field; condition, mask, dims)

struct Integral end

function Reduction(int::Integral, field::AbstractField; condition = nothing, mask = 0, dims)
    dims = dims isa Colon ? (1, 2, 3) : tupleit(dims)
    dx = reduction_grid_metric(dims)
    return Reduction(sum!, condition_operand(field * dx, condition, mask), dims)
end

"""
    Integral(field; dims=:)

Return a `Reduction` representing a spatial integral of `field` over `dims`.
"""
Integral(field::AbstractField; condition = nothing, mask = 0, dims=:) = Reduction(Integral(), condition_operand(field, condition, mask), dims)

#####
##### show
#####

Base.summary(r::Reduction{<:Average}) = string("Average of ", summary(r.operand), " over dims ", r.dims)
Base.summary(r::Reduction{<:Integral}) = string("Integral of ", summary(r.operand), " over dims ", r.dims)


text:
"""
MkCell provides functions for calculating optimal supercells.

Based on the description of the approach given by
    > B. Militzer, High Energy Density Phys. 21 (2016) 8–15. doi:10.1016/j.hedp.2016.09.003
"""
module MkCell

export cellopt, cell_abc, cell_angles

using LinearAlgebra, StaticArrays

include("vecutils.jl")
include("scoring.jl")

"""
Find an optimal supercell by scaling `cell` by `m` unit cells, searching over
a 3x3 supercell matrix where each element ranges from `-n` to `+n`.

This function follows the approach described by
    > B. Militzer, High Energy Density Phys. 21 (2016) 8–15. doi:10.1016/j.hedp.2016.09.003

An optimal cell is chosen according to the following four criteria:
1. Minimizing the sphere radius required to contain the supercell
2. Maximizing the periodic distance between images
3. Minimizing the distortion of the supercell angles from 90 degrees
4. Minimizing the difference in lattice vector length
"""
function cellopt(cell, m, n; verbose=false)
    Vsuper = abs(det(cell)*m)  # target supercell volume
    a, b, c = SVector{3}(cell[1,:]), SVector{3}(cell[2,:]), SVector{3}(cell[3,:])
    count = 0;
    best_cell = MMatrix{3,3}(Array{Float64}(undef, 3, 3))
    # The first value should be set to infinity so that it will always be
    # higher than the current cell for the first iteration.
    best_cell_score = [Inf NaN NaN NaN]
    cur_cell = MMatrix{3,3}(Array{Float64}(undef, 3, 3));
    best_matrix = MMatrix{3,3}(Array{Float64}(undef, 3, 3));
    cur_matrix = MMatrix{3,3}(Array{Float64}(undef, 3, 3));
    asuper = MVector{3}(Array{Float64}(undef, 3));
    bsuper = MVector{3}(Array{Float64}(undef, 3));
    csuper = MVector{3}(Array{Float64}(undef, 3));
    for ia=UnitRange{Float64}(1,n), ja=UnitRange{Float64}(-n,n), ka=UnitRange{Float64}(-n,n)
        asuper[:] = ia*a + ja*b + ka*c;
        if norm(asuper) <= 0
            continue
        end # if
        for ib=UnitRange{Float64}(-ia,ia), jb=UnitRange{Float64}(-n,n), kb=UnitRange{Float64}(-n,n)
            bsuper[:] = ib*a + jb*b + kb*c;
            if norm(bsuper) <= 0
                continue
            end # if

            aib = abs(ib)
            ajb = abs(jb)
            for ic=UnitRange{Float64}(-aib,aib), jc=UnitRange{Float64}(-ajb,ajb), kc=UnitRange{Float64}(-n,n)
                csuper[:] = ic*a + jc*b + kc*c;
                cur_matrix[1] = ia; cur_matrix[2] = ib; cur_matrix[3] = ic;
                cur_matrix[4] = ja; cur_matrix[5] = jb; cur_matrix[6] = jc;
                cur_matrix[7] = ka; cur_matrix[8] = kb; cur_matrix[9] = kc;
                cur_cell[:,:] = [asuper'; bsuper'; csuper']
                Vsuper_calc = det(cur_cell);
                if abs(Vsuper_calc) != Vsuper
                    continue
                end # if
                # This is a valid supercell, now try to determine if it's the optimal cell so far
                cellchoose!(best_cell, best_cell_score, cur_cell, best_matrix, cur_matrix)
                count = count + 1;
            end # for
        end # for
    end # for
    # reverse the sign of dmin so that it's positive, consistent with reality
    best_cell_score[2] = -best_cell_score[2]
    if verbose
        as = best_cell[1,:]
        bs = best_cell[2,:]
        cs = best_cell[3,:]
        println("Best cell", best_cell)
        println("a b c ", norm(as), " ", norm(bs), " ", norm(cs))
        println("α β γ (°) ", ang(as, bs)*180/π, " ", ang(as, cs)*180/π, " ", ang(bs, cs)*180/π)
        println("Matrix ", best_matrix)
        println("Rmax ", best_cell_score[1])
        println("dmin ", best_cell_score[2])
        println("Angle deviation ", best_cell_score[3])
        println("Length deviation ", best_cell_score[4])
        println("Cells considered ", count)
    end
    return best_cell, best_matrix, best_cell_score, count
end # function

# TODO: add a function that can calculate optimal cells over a range of unit cells
# multiples to determine which multiples have dmin larger than all smaller cells

# TODO: benchmark for maximum performance.
end # module


text:
# Wisdom & Hernandez version of Kepler solver, but with quartic convergence.

function cubic1(a::T, b::T, c::T) where {T <: Real}
a3 = a*third
Q = a3^2 - b*third
R = a3^3 + 0.5*(-a3*b + c)
if R^2 < Q^3
# println("Error in cubic solver ",R^2," ",Q^3)
 return -c/b
else
 A = -sign(R)*cbrt(abs(R) + sqrt(R*R - Q*Q*Q))
 if A == 0.0
   B = 0.0
  else
   B = Q/A
 end
 x1 = A + B - a3
 return x1
end
return
end

function calc_ds_opt(y::T,yp::T,ypp::T,yppp::T) where {T <: Real}
# Computes quartic Newton's update to equation y=0 using first through 3rd derivatives.
# Uses techniques outlined in Murray & Dermott for Kepler solver.
# Rearrange to reduce number of divisions:
num = y*yp
den1 = yp*yp-y*ypp*.5
den12 = den1*den1
den2 = yp*den12-num*.5*(ypp*den1-third*num*yppp)
return -y*den12/den2
end

function solve_kepler!(h::T,k::T,x0::Array{T,1},v0::Array{T,1},beta0::T,
   r0::T,s0::T,state::Array{T,1}) where {T <: Real}
# Solves elliptic Kepler's equation for both elliptic and hyperbolic cases.
# Initial guess (if s0 = 0):
r0inv = inv(r0)
beta0inv = inv(beta0)
signb = sign(beta0)
sqb = sqrt(signb*beta0)
zeta = k-r0*beta0
#eta = dot(x0,v0)
eta = x0[1]*v0[1]+x0[2]*v0[2]+x0[3]*v0[3]
sguess = 0.0
if s0 == zero(T)
  # Use cubic estimate:
  if zeta != zero(T)
    sguess = cubic1(3eta/zeta,6r0/zeta,-6h/zeta)
  else
    if eta != zero(T)
      reta = r0/eta
      disc = reta^2+2h/eta
      if disc > zero(T)
        sguess =-reta+sqrt(disc)
      else
       sguess = h*r0inv
      end
    else
      sguess = h*r0inv
    end
  end
  s = copy(sguess)
else
  s = copy(s0)
end
s0 = copy(s)
y = zero(T); yp = one(T)
iter = 0
ds = Inf
KEPLER_TOL = sqrt(eps(h))
ITMAX = 20
while iter == 0 || (abs(ds) > KEPLER_TOL && iter < ITMAX)
  xx = sqb*s
  if beta0 > 0
    #sx = sin(xx); cx = cos(xx)
    sx,cx = sincos(xx)
  else
    cx = cosh(xx); sx = exp(xx)-cx
  end
  sx *= sqb
# Third derivative:
  yppp = zeta*cx - signb*eta*sx
# First derivative:
  yp = (-yppp+ k)*beta0inv
# Second derivative:
  ypp = signb*zeta*beta0inv*sx + eta*cx
  y  = (-ypp + eta +k*s)*beta0inv - h  # eqn 35
# Now, compute fourth-order estimate:
#  ds = calc_ds_opt(y,yp,ypp,yppp)
  ds = -y/yp
  s += ds
  iter +=1
end
if iter == ITMAX
  println("Reached max iterations in solve_kepler: h ",h," s0: ",s0," s: ",s," ds: ",ds)
end
#println("sguess: ",sguess," s: ",s," s-sguess: ",s-sguess," ds: ",ds," iter: ",iter)
# Since we updated s, need to recompute:
xx = 0.5*sqb*s
if beta0 > 0
  #sx = sin(xx); cx = cos(xx)
  sx,cx = sincos(xx)
else
  cx = cosh(xx); sx = exp(xx)-cx
end
# Now, compute final values:
g1bs = 2.0*sx*cx/sqb
g2bs = 2.0*signb*sx^2*beta0inv
f = one(T) - k*r0inv*g2bs # eqn (25)
g = r0*g1bs + eta*g2bs # eqn (27)
for j=1:3
# Position is components 2-4 of state:
  state[1+j] = x0[j]*f+v0[j]*g
end
r = sqrt(state[2]*state[2]+state[3]*state[3]+state[4]*state[4])
rinv = inv(r)
dfdt = -k*g1bs*rinv*r0inv
dgdt = (r0-r0*beta0*g2bs+eta*g1bs)*rinv
for j=1:3
# Velocity is components 5-7 of state:
  state[4+j] = x0[j]*dfdt+v0[j]*dgdt
end
return s,f,g,dfdt,dgdt,cx,sx,g1bs,g2bs,r,rinv,ds,iter
end

function kep_ell_hyp!(x0::Array{T,1},v0::Array{T,1},r0::T,k::T,h::T,
  beta0::T,s0::T,state::Array{T,1}) where {T <: Real}
# Solves equation (35) from Wisdom & Hernandez for the elliptic case.
# Now, solve for s in elliptical Kepler case:
f = zero(T); g=zero(T); dfdt=zero(T); dgdt=zero(T); cx=zero(T);sx=zero(T);g1bs=zero(T);g2bs=zero(T)
s=zero(T); ds = zero(T); r = zero(T);rinv=zero(T); iter=0
if beta0 > zero(T) || beta0 < zero(T)
   s,f,g,dfdt,dgdt,cx,sx,g1bs,g2bs,r,rinv,ds,iter = solve_kepler!(h,k,x0,v0,beta0,r0,
    s0,state)
else
  println("Not elliptic or hyperbolic ",beta0," x0 ",x0)
  r= zero(T); fill!(state,zero(T)); rinv=zero(T); s=zero(T); ds=zero(T); iter = 0
end
state[8]= r
state[9] = (state[2]*state[5]+state[3]*state[6]+state[4]*state[7])*rinv
# recompute beta:
# beta is element 10 of state:
state[10] = 2.0*k*rinv-(state[5]*state[5]+state[6]*state[6]+state[7]*state[7])
# s is element 11 of state:
state[11] = s
# ds is element 12 of state:
state[12] = ds
return iter
end

function kep_ell_hyp!(x0::Array{T,1},v0::Array{T,1},r0::T,k::T,h::T,
  beta0::T,s0::T,state::Array{T,1},jacobian::Array{T,2}) where {T <: Real}
# Computes the Jacobian as well
# Solves equation (35) from Wisdom & Hernandez for the elliptic case.
r0inv = inv(r0)
beta0inv = inv(beta0)
# Now, solve for s in elliptical Kepler case:
f = zero(T); g=zero(T); dfdt=zero(T); dgdt=zero(T); cx=zero(T);sx=zero(T);g1bs=zero(T);g2bs=zero(T)
s=zero(T); ds=zero(T); r = zero(T);rinv=zero(T); iter=0
if beta0 > zero(T) || beta0 < zero(T)
   s,f,g,dfdt,dgdt,cx,sx,g1bs,g2bs,r,rinv,ds,iter = solve_kepler!(h,k,x0,v0,beta0,r0,
    s0,state)
# Compute the Jacobian.  jacobian[i,j] is derivative of final state variable q[i]
# with respect to initial state variable q0[j], where q = {x,v} & q0 = {x0,v0}.
  fill!(jacobian,zero(T))
  compute_jacobian!(h,k,x0,v0,beta0,s,f,g,dfdt,dgdt,cx,sx,g1bs,g2bs,r0,r,jacobian)
else
  println("Not elliptic or hyperbolic ",beta0," x0 ",x0)
  r= zero(T); fill!(state,zero(T)); rinv=zero(T); s=zero(T); ds=zero(T); iter = 0
end
# recompute beta:
state[8]= r
state[9] = (state[2]*state[5]+state[3]*state[6]+state[4]*state[7])*rinv
# beta is element 10 of state:
state[10] = 2.0*k*rinv-(state[5]*state[5]+state[6]*state[6]+state[7]*state[7])
# s is element 11 of state:
state[11] = s
# ds is element 12 of state:
state[12] = ds
return iter
end

function compute_jacobian!(h::T,k::T,x0::Array{T,1},v0::Array{T,1},beta0::T,s::T,
  f::T,g::T,dfdt::T,dgdt::T,cx::T,sx::T,g1::T,g2::T,r0::T,r::T,jacobian::Array{T,2}) where {T <: Real}
# Compute the Jacobian.  jacobian[i,j] is derivative of final state variable q[i]
# with respect to initial state variable q0[j], where q = {x,v,k} & q0 = {x0,v0,k}.
# Now, compute the Jacobian: (9/18/2017 notes)
g0 = one(T)-beta0*g2
g3 = (s-g1)/beta0
eta = dot(x0,v0)  # unnecessary to divide by r0 for multiply for \dot\alpha_0
absv0 = sqrt(dot(v0,v0))
dsdbeta = (2h-r0*(s*g0+g1)+k/beta0*(s*g0-g1)-eta*s*g1)/(2beta0*r)
dsdr0 = -(2k/r0^2*dsdbeta+g1/r)
dsda0 = -g2/r
dsdv0 = -2absv0*dsdbeta
dsdk = 2/r0*dsdbeta-g3/r
dbetadr0 = -2k/r0^2
dbetadv0 = -2absv0
dbetadk  = 2/r0
# "p" for partial derivative:
for i=1:3
  pxpr0[i] = k/r0^2*g2*x0[i]+g1*v0[i]
  pxpa0[i] = g2*v0[i]
  pxpk[i]  = -g2/r0*x0[i]
  pxps[i]  = -k/r0*g1*x0[i]+(r0*g0+eta*g1)*v0[i]
  pxpbeta[i] = -k/(2beta0*r0)*(s*g1-2g2)*x0[i]+1/(2beta0)*(s*r0*g0-r0*g1+s*eta*g1-2*eta*g2)*v0[i]
  prvpr0[i] = k*g1/r0^2*x0[i]+g0*v0[i]
  prvpa0[i] = g1*v0[i]
  prvpk[i] = -g1/r0*x0[i]
  prvps[i] = -k*g0/r0*x0[i]+(-beta0*r0*g1+eta*g0)*v0[i]
  prvpbeta[i] = -k/(2beta0*r0)*(s*g0-g1)*x0[i]+1/(2beta0)*(-s*r0*beta0*g1+eta*s*g0-eta*g1)*v0[i]
end
prpr0 = g0
prpa0 = g1
prpk  = g2
prps = (k-beta0*r0)*g1+eta*g0
prpbeta = 1/(2beta0)*(s*(k-beta0*r0)*g1+eta*s*g0-eta*g1-2k*g2)
for i=1:3
  dxdr0[i] = pxps[i]*dsdr0 + pxpbeta[i]*dbetadr0 + pxpr0[i]
  dxda0[i] = pxps[i]*dsda0 + pxpa0[i]
  dxdv0[i] = pxps[i]*dsdv0 + pxpbeta[i]*dbetadv0
  dxdk[i]  = pxps[i]*dsdk  + pxpbeta[i]*dbetadk + pxpk[i]
  drvdr0[i] = prvps[i]*dsdr0 + prvpbeta[i]*dbetadr0 + prvpr0[i]
  drvda0[i] = prvps[i]*dsda0 + prvpa0[i]
  drvdv0[i] = prvps[i]*dsdv0 + prvpbeta[i]*dbetadv0
  drvdk[i]  = prvps[i]*dsdk  + prvpbeta[i]*dbetadk +prvpk[i]
end
drdr0 = prpr0 + prps*dsdr0 + prpbeta*dbetadr0
drda0 = prpa0 + prps*dsda0
drdv0 = prps*dsdv0 + prpbeta*dbetadv0
drdk  = prpk + prps*dsdk + prpbeta*dbetadk
for i=1:3
  vtmp[i] = dfdt*x0[i]+dgdt*v0[i]
  dvdr0[i] = (drvdr0[i]-drdr0*vtmp[i])/r
  dvda0[i] = (drvda0[i]-drda0*vtmp[i])/r
  dvdv0[i] = (drvdv0[i]-drdv0*vtmp[i])/r
  dvdk[i]  = (drvdk[i] -drdk *vtmp[i])/r
end
# Now, compute Jacobian:
for i=1:3
  jacobian[  i,  i] = f
  jacobian[  i,3+i] = g
  jacobian[3+i,  i] = dfdt
  jacobian[3+i,3+i] = dgdt
  jacobian[  i,7] = dxdk[i]
  jacobian[3+i,7] = dvdk[i]
end
for j=1:3
  for i=1:3
    jacobian[  i,  j] += dxdr0[i]*x0[j]/r0 + dxda0[i]*v0[j]
    jacobian[  i,3+j] += dxdv0[i]*v0[j]/absv0 + dxda0[i]*x0[j]
    jacobian[3+i,  j] += dvdr0[i]*x0[j]/r0 + dvda0[i]*v0[j]
    jacobian[3+i,3+j] += dvdv0[i]*v0[j]/absv0 + dvda0[i]*x0[j]
  end
end
jacobian[7,7]=one(T)
return
end


text:
export StateTransformedEnv

struct StateTransformedEnv{P,M,E<:AbstractEnv} <: AbstractEnvWrapper
    env::E
    state_mapping::P
    state_space_mapping::M
end

"""
    StateTransformedEnv(env; state_mapping=identity, state_space_mapping=identity)

`state_mapping` will be applied on the original state when calling `state(env)`,
and similarly `state_space_mapping` will be applied when calling `state_space(env)`.
"""
StateTransformedEnv(env; state_mapping=identity, state_space_mapping=identity) =
    StateTransformedEnv(env, state_mapping, state_space_mapping)

RLBase.state(env::StateTransformedEnv, args...; kwargs...) =
    env.state_mapping(state(env.env, args...; kwargs...))

RLBase.state_space(env::StateTransformedEnv, args...; kwargs...) =
    env.state_space_mapping(state_space(env.env, args...; kwargs...))


text:
using ProgressBars, Plots, StatsBase, Statistics, LaTeXStrings, JLD

function resultsPart1(N)
    return rand(0:9,N)
end

function resultsPart2(N)
    TotRandoms = rand(0:9,N)
    TotRandNumbers = [length(findall(x->x==n,TotRandoms)) for n ∈ 0:9]
    return std(TotRandNumbers)/N
end

function Linear_fit(X, Y)
    A = [hcat(X) reshape(ones(length(X)), length(X), 1)]
    b = reshape(hcat(Y), length(Y), 1)
    line = (A \ b)
    return line
end

N = 10^6
TotRandoms = resultsPart1(N)
save("../../Data/Q4/Q4-hist.jld", "Data", TotRandoms)

plot([-0.5,9.5],[N/10,N/10], linestyle = :dash, label = L"Y = \frac{N}{10}", c = :black, linewidth = 4, ylims = (-5000,1.35*10^5))
histogram!(TotRandoms, bins = -0.5:9.5, xticks = 0:9, label = L"Histogram", framestyle = :box, c = :steelblue)
plot!(xlabel = L"Number", ylabel = L"Count",title = L"%$N\ Random\ Numbers\ in\ [0, 9]")
savefig("../../Figs/Q4/Q4-Hist.pdf")

#=====================##=====================##=====================#

span_ = ceil.(Int,exp.(2:0.25:10))
Data = zeros(1000,length(span_))

for n in ProgressBar(1:1000)
    Data[n,:] = [resultsPart2(N) for N in span_]
end

CVList = [mean(Data[:,n]) for n ∈ 1:length(span_)]
Line = Linear_fit(log.(span_),log.(CVList))
X = log.(span_)
Y = X .* Line[1,1] .+ Line[2,1]

save("../../Data/Q4/Q4-CV.jld", "Data", CVList)

plot(X,Y, label = L"Y = %$(round(Line[1],digits=3))\ X + %$(round(Line[2],digits=3))", line = :dash, c = :black)
scatter!(log.(span_), log.(CVList),label = L"Data\ Point", framestyle = :box, c = :steelblue)
scatter!(title = L"log-log\ Plot\ CV-N\ (e^2>N>e^{10},\ 10000\ runs)",
    ylabel = L"ln_{(CV_{(N)})}", xlabel = L"ln_{(N)}")
savefig("../../Figs/Q4/Q4-Scat.pdf")


text:
function corner_coordinates(tess::Tessellation)
    # Each cell needs the number of points + line from last to first point + NaN
    n_points = mapreduce(x -> length(x) + 2, +, tess.Cells)
    p = Vector{GeometryBasics.Point2{Float64}}(undef, n_points)

    count = 0
    for cell in tess.Cells
        for corner in cell
            count += 1
            p[count] = GeometryBasics.Point2(getx(corner), gety(corner))
        end

        count += 1
        p[count] = GeometryBasics.Point2(getx(cell[1]), gety(cell[1]))

        count += 1
        p[count] = GeometryBasics.Point2(NaN, NaN)
    end

    return p
end

@recipe function plot(tess::Tessellation)
    @series begin
        xlims --> (left(tess.EnclosingRectangle), right(tess.EnclosingRectangle))
        ylims --> (lower(tess.EnclosingRectangle), upper(tess.EnclosingRectangle))
        seriestype --> :path
        label --> "Voronoi cells"

        corner_coordinates(tess)
    end
end


text:
using Oceananigans.Grids: xnode, znode
using Oceananigans.TimeSteppers: update_state!

function run_rayleigh_benard_regression_test(arch, grid_type)

    #####
    ##### Parameters
    #####
          α = 2                 # aspect ratio
          n = 1                 # resolution multiple
         Ra = 1e6               # Rayleigh number
    Nx = Ny = 8n * α            # horizontal resolution
    Lx = Ly = 1.0 * α           # horizontal extent
         Nz = 16n               # vertical resolution
         Lz = 1.0               # vertical extent
         Pr = 0.7               # Prandtl number
          a = 1e-1              # noise amplitude for initial condition
         Δb = 1.0               # buoyancy differential

    # Rayleigh and Prandtl determine transport coefficients
    ν = sqrt(Δb * Pr * Lz^3 / Ra)
    κ = ν / Pr

    #####
    ##### Model setup
    #####

    if grid_type == :regular
        grid = RegularRectilinearGrid(size=(Nx, Ny, Nz), extent=(Lx, Ly, Lz))
    elseif grid_type == :vertically_unstretched
        zF = range(-Lz, 0, length=Nz+1)
        grid = VerticallyStretchedRectilinearGrid(architecture=arch, size=(Nx, Ny, Nz), x=(0, Lx), y=(0, Ly), z_faces=zF)
    end

    # Force salinity as a passive tracer (βS=0)
    c★(x, z) = exp(4z) * sin(2π/Lx * x)
    Fc(i, j, k, grid, clock, model_fields) = 1/10 * (c★(xnode(Center(), i, grid), znode(Center(), k, grid)) - model_fields.c[i, j, k])

    bbcs = FieldBoundaryConditions(top = BoundaryCondition(Value, 0.0),
                                   bottom = BoundaryCondition(Value, Δb))

    model = IncompressibleModel(
               architecture = arch,
                       grid = grid,
                    closure = IsotropicDiffusivity(ν=ν, κ=κ),
                    tracers = (:b, :c),
                   buoyancy = Buoyancy(model=BuoyancyTracer()),
        boundary_conditions = (b=bbcs,),
                    forcing = (c=Forcing(Fc, discrete_form=true),)
    )

    # Lz/Nz will work for both the :regular and :vertically_unstretched grids.
    Δt = 0.01 * min(model.grid.Δx, model.grid.Δy, Lz/Nz)^2 / ν

    # We will manually change the stop_iteration as needed.
    simulation = Simulation(model, Δt=Δt, stop_iteration=0)

    # The type of the underlying data, not the offset array.
    ArrayType = typeof(model.velocities.u.data.parent)

    spinup_steps = 1000
      test_steps = 100

    prefix = "rayleigh_benard"

    checkpointer = Checkpointer(model, schedule=IterationInterval(test_steps), prefix=prefix,
                                dir=joinpath(dirname(@__FILE__), "data"))

    #####
    ##### Initial condition and spinup steps for creating regression test data
    #####

    #=
    @warn "Generating new data for the Rayleigh-Benard regression test."

    ξ(z) = a * rand() * z * (Lz + z) # noise, damped at the walls
    b₀(x, y, z) = (ξ(z) - z) / Lz
    set!(model, b=b₀)

    simulation.stop_iteration = spinup_steps-test_steps
    run!(simulation)

    push!(simulation.output_writers, checkpointer)
    simulation.stop_iteration += 2test_steps
    run!(simulation)
    =#

    #####
    ##### Regression test
    #####

    # Load initial state
    initial_filename = joinpath(dirname(@__FILE__), "data", prefix * "_iteration$spinup_steps.jld2")

    solution₀, Gⁿ₀, G⁻₀ = get_fields_from_checkpoint(initial_filename)

    model.velocities.u.data.parent .= ArrayType(solution₀.u)
    model.velocities.v.data.parent .= ArrayType(solution₀.v)
    model.velocities.w.data.parent .= ArrayType(solution₀.w)
    model.tracers.b.data.parent    .= ArrayType(solution₀.b)
    model.tracers.c.data.parent    .= ArrayType(solution₀.c)

    model.timestepper.Gⁿ.u.data.parent .= ArrayType(Gⁿ₀.u)
    model.timestepper.Gⁿ.v.data.parent .= ArrayType(Gⁿ₀.v)
    model.timestepper.Gⁿ.w.data.parent .= ArrayType(Gⁿ₀.w)
    model.timestepper.Gⁿ.b.data.parent .= ArrayType(Gⁿ₀.b)
    model.timestepper.Gⁿ.c.data.parent .= ArrayType(Gⁿ₀.c)

    model.timestepper.G⁻.u.data.parent .= ArrayType(G⁻₀.u)
    model.timestepper.G⁻.v.data.parent .= ArrayType(G⁻₀.v)
    model.timestepper.G⁻.w.data.parent .= ArrayType(G⁻₀.w)
    model.timestepper.G⁻.b.data.parent .= ArrayType(G⁻₀.b)
    model.timestepper.G⁻.c.data.parent .= ArrayType(G⁻₀.c)

    model.clock.iteration = spinup_steps
    model.clock.time = spinup_steps * Δt
    length(simulation.output_writers) > 0 && pop!(simulation.output_writers)

    # Step the model forward and perform the regression test
    update_state!(model)

    for n in 1:test_steps
        time_step!(model, Δt, euler=false)
    end

    final_filename = joinpath(@__DIR__, "data", prefix * "_iteration$(spinup_steps+test_steps).jld2")

    solution₁, Gⁿ₁, G⁻₁ = get_fields_from_checkpoint(final_filename)

    test_fields =  CUDA.@allowscalar (u = Array(interior(model.velocities.u)),
                                      v = Array(interior(model.velocities.v)),
                                      w = Array(interior(model.velocities.w)[:, :, 1:Nz]),
                                      b = Array(interior(model.tracers.b)),
                                      c = Array(interior(model.tracers.c)))

    correct_fields = (u = Array(interior(solution₁.u, model.grid)),
                      v = Array(interior(solution₁.v, model.grid)),
                      w = Array(interior(solution₁.w, model.grid)),
                      b = Array(interior(solution₁.b, model.grid)),
                      c = Array(interior(solution₁.c, model.grid)))

    summarize_regression_test(test_fields, correct_fields)

    CUDA.allowscalar(true)
    @test all(test_fields.u .≈ correct_fields.u)
    @test all(test_fields.v .≈ correct_fields.v)
    @test all(test_fields.w .≈ correct_fields.w)
    @test all(test_fields.b .≈ correct_fields.b)
    @test all(test_fields.c .≈ correct_fields.c)
    CUDA.allowscalar(false)

    return nothing
end


text:
__precompile__(false)

module CameraManager

using Distributed
using GLFW
using ModernGL

using MathManager
using TimeManager
using WindowManager

function reset()
  global lastCursorPos = [0,0]
  global cursorPos_old = [0,0]
  global cursorPos = [0,0]
  global mouseMove = [0,0]
end

reset()

#MVP = eye(Mat4x4f)
#MVP *= transform([0.0f0,0.0f0,0.0f0],[0.0f0,0.0f0,0.0f0],Float32[])

keyFB = 0
keyLR = 0
keyUD = 0
speed = false

mouseKeyPressed = false
keyPressed = false
keyPressing = false
keyValue = 0

getKey() = keyValue, keyPressed

export getKey

resetKeys() = global keyPressed = false

export resetKeys

function cleanUp()
  reset()
  resetKeys()
end

"""
camera object with holds position, rotation, scaling and various matrices like MVP
"""
mutable struct Camera
  moved::Bool

  position::AbstractArray
  rotation::AbstractArray
  scale::AbstractArray

  viewMat::AbstractArray
  projectionMat::AbstractArray
  modelMat::AbstractArray

  translateMat::AbstractArray
  rotationMat::AbstractArray
  scalingMat::AbstractArray

  MVP::AbstractArray

  Camera() = new(true, zerosVector3f,zerosVector3f,onesVector3f,
  eyeMat4x4f,eyeMat4x4f,eyeMat4x4f,eyeMat4x4f,eyeMat4x4f,eyeMat4x4f,eyeMat4x4f)
end

export Camera

CAMERA = Camera()

"""
gets forward vector of camera direction
"""
forward(camera::Camera) = MathManager.forward(camera.rotationMat)

export forward

"""
gets right vector of camera direction
"""
right(camera::Camera) = MathManager.right(camera.rotationMat)

export right

"""
gets up vector of camera direction
"""
up(camera::Camera) = MathManager.up(camera.rotationMat)

export up

"""
sets projection matrix
"""
setProjection(camera::Camera, m::AbstractArray) = (camera.projectionMat = m)

export setProjection

"""
sets view matrix
"""
setView(camera::Camera, m::AbstractArray) = (camera.viewMat = m)

export setView

VIEW_KEYS=false

"""
event which catches keyboard inputs.
here keys for wireframe, fullscreen and camera movement are defined
"""
function OnKey(window::Window, key::Number, scancode::Number, action::Number, mods::Number, unicode::Number)
  #println("K")
  #key, scancode, action, mods, unicode = Int32(key), Int32(scancode), Int32(action), Int32(mods), Int32(unicode)

  if key == 70 && action == 1 # f
    # ...
  elseif key == 71 && action == 1 # g
    #global wireframe
    #wireframe = !wireframe
    #println("Wireframe: ",wireframe)

  elseif key == 72 && action == 1 # h
    #global fullscreen=!fullscreen

  elseif key == 87 #w
    global keyFB = (action > 0) ? 1 : 0
  elseif key == 83 #s
    global keyFB = (action > 0) ? -1 : 0
  elseif key == 65 #a
    global keyLR = (action > 0) ? -1 : 0
  elseif key == 68 #d
    global keyLR = (action > 0) ? 1 : 0
  elseif key == 32 #space
    global keyUD = (action > 0) ? 1 : 0
  elseif key == 67 || key == 341 # c || lctrl
    global keyUD = (action > 0) ? -1 : 0
  elseif key == 340 #lshift
    global speed = (action > 0)
  end

  if key == 75 && action == 1 #k
    global VIEW_KEYS = !VIEW_KEYS
    println("VIEW_KEYS = $VIEW_KEYS")
  end

  if VIEW_KEYS println("K: ", key, " A: ", action, " C: ", scancode, " M: ", mods) end

  if action == 1 global keyPressed = true end
  global keyPressing = action != 0
  global keyValue = key
  nothing
end

export OnKey

"""
event which catches mouse key inpits and hides/shows cursor when mouse button is pressed
"""
function OnMouseKey(window::Window, key::Number, action::Number, mods::Number)
  #println("M")

  key = Int32(key)
  action = Int32(action)
  mods = Int32(mods)

  if keyPressing return end

  global cursorPos_old, cursorPos
  global mouseKeyPressed = action != 0

  if action == 1
    cursorPos_old = cursorPos
    GLFW.SetInputMode(window.ref, GLFW.CURSOR, GLFW.CURSOR_DISABLED)
  end
  if action == 0
    GLFW.SetInputMode(window.ref, GLFW.CURSOR, GLFW.CURSOR_NORMAL)
  end
  nothing
end

export OnMouseKey

"""
event which catches mouse position for camera rotation event
"""
function OnCursorPos(window::Window, pos::Tuple{Number,Number}) #window, x::Number, y::Number)
  global CAMERA, lastCursorPos, cursorPos, mouseMove, mouseKeyPressed
  #println("C")
  (x,y) = pos
  t = [x/800f0,y/600f0]
  cursorPos = t
  l = lastCursorPos; if l == [] l = t end
  n = l - t
  lastCursorPos = t
  mouseMove = n

  if mouseKeyPressed OnRotate(CAMERA) end
  nothing
end

export OnCursorPos

"""
rotates camera
"""
function rotate(camera::Camera, rotation::AbstractArray)
  camera.rotation += rotation
  camera.rotationMat = computeRotation(camera.rotation)
  camera.moved=true
  #if length(procs()) > 1 put!(channels[:CAMERA],(camera.position,camera.rotation)) end
  #println("rotated")
end

export rotate

"""
moves camera, adds vector to current position
"""
function move(camera::Camera, position::AbstractArray)
  camera.position += position
  camera.translateMat = translation(camera.position)
  camera.moved=true
  #if length(procs()) > 1 put!(channels[:CAMERA],(camera.position,camera.rotation)) end
  #println("moved")
end

export move

"""
event which calculates cursor position shifts and calls rotate function
"""
function OnRotate(camera::Camera)
  global cursorPos, cursorPos_old
  mx = cursorPos[1] - cursorPos_old[1]
  my = cursorPos[2] - cursorPos_old[2]
  #println(mx)
  cursorPos_old = cursorPos
  rotate(camera, [-mx*2f0,my*2f0,0f0]) #[-mx*2,my*2,0f0] #Vec3f((-mx+0.5f0),(-0.5f0+my),0f0)
end

export OnRotate

oldposition = CAMERA.position
shiftposition = [0.0f0,0,0]

"""
sets camera position
"""
function setPosition(camera::Camera, position::AbstractArray)
  camera.position = position
  global oldposition = position
  camera.translateMat = translation(camera.position)
  camera.moved = true
  if length(procs()) > 1 put!(channels[:CAMERA],(camera.position,camera.rotation)) end
end

export setPosition

"""
event which updates positions shifts (left,right,up,down,forward,back)
key is (left,right,up,down,forward,back)
m is direction value with weight (positive, negative)
"""
function OnMove(camera::Camera, key::Symbol, f::Number)
  global shiftposition
  if key == :FORWARD  move(camera, forward(camera)*f)
  elseif key == :RIGHT  move(camera, right(camera)*-f) #+Vec3f(-right*0.02f0,-up*0.02f0,forward*0.02f0)
  elseif key == :UP  move(camera, up(camera)*-f)
  end

  #dif = camera.position
  #step = trunc.(dif)
  #shiftposition = dif - step
  #step = abs.(step)

  #if step[1] >= 2 camera.position[1] = 0 end
  #if step[2] >= 2 camera.position[2] = 0 end
  #if step[3] >= 2 camera.position[3] = 0 end

  #  dist = sqrt(camera.position[1]^2+camera.position[2]^2+camera.position[3]^2)
  # if dist >= 1 camera.position = [0f0,0,0] end
end

export OnMove

"""
update function where camera translation is update only when camera was moved by input.
here cameras MVP Matrix is updated aswell
"""
function Update(camera::Camera)
  camera.viewMat = camera.scalingMat * camera.rotationMat * camera.translateMat
  camera.MVP = camera.modelMat * camera.projectionMat * camera.viewMat
end

export Update

"""
event which is called by game loop and calls real update function
this event resets camera moved state
"""
function OnUpdate(camera::Camera)
  #if isFocus return end

  if OnTime(0.01)
    f=GetTimePassed()*(!speed ? 5 : 100)
    if keyFB != 0 OnMove(camera, :FORWARD, keyFB*f) end
    if keyLR != 0 OnMove(camera, :RIGHT, keyLR*f) end
    if keyUD != 0 OnMove(camera, :UP, keyUD*f) end
  end

  if camera.moved
    Update(camera)
  end

  r = camera.moved
  camera.moved = false
  r
end

end #CameraManager

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
"""
    KnetModel

    * tensors: A dictionary. Given the tensor name as a string, returns the actual tensor.

    * model_layers: returns the list of layers (the actual layers themselves)

    * model_inputs: returns the list of inputs (the names of the tensors)

    * model_outputs: returns the list of outputs (the names of the outputs)
"""

mutable struct KnetModel
    tensors #dictionary: string -> tensor
    model_layers
    model_inputs
    model_outputs
    graph
end

function KnetModel(g::KnetOnnx.Types.Graph)
    """
    KnetModel constructor.
    Uses the Graph.
    """
    model_layers = get_ModelLayers(g)
    tensors = TensorDict(model_layers, g)
    model_inputs = [i.name for i in g.input]
    model_outputs = [o.name for o in g.output]
    graph = g
    KnetModel(tensors, model_layers, model_inputs, model_outputs, graph)
end

function KnetModel(path::String)
    """
    KnetModel constructor #2
    Uses the path of the ONNX file.
    """
    KnetModel(ONNXtoGraph(path))
end

"""
    TensorDict
    Initializes KnetModel.tensors by putting the names in as keys, but the values are Nothing.
"""
#omitted weights, might need them later
function TensorDict(model_layers, g)
    tensors = Dict()
    for layer in model_layers
        for input in layer.inputs; tensors[input] = Nothing; end
        for input in layer.outputs; tensors[input] = Nothing; end
    end
    for t in keys(g.initializer)
        if t in keys(tensors); tensors[t] = g.initializer[t]; end
    end
    tensors
end

"""
    ModelLayer
        * inputs: a list of tensor names.
          These tensors will be used for forward calculation of the layer.
        * outputs: a list of tensor names.
          The outputs of the forward calculation will be saved to Model.tensors under these keys.
        * layer: a Knet Layer.
          If you are constructing your own ModelLayer make sure the number of inputs and outputs matches the functionality of the KnetLayer you are using.

"""
mutable struct ModelLayer
    inputs #list of strings
    layer # a KnetLayer
    outputs #list of strings
end

function ModelLayer(node, g)
    (args, layer, outputs) = convert(node, g)
    ModelLayer(args, layer, outputs)
end

function get_ModelLayers(g)
    ModelLayers = []
    for node in g.node; push!(ModelLayers, ModelLayer(node, g)); end
    return ModelLayers
end

# FORWARD CALCULATIONS
function forward(km::KnetModel, ml::ModelLayer)

        # GATHER INPUTS
    for input in ml.inputs
        if km.tensors[input] == Nothing; return "oops!"; end
    end

        # FORWARD PASS
        # if only one input is requried, pass the first element
        # if more than one input is required, pass all elements
        # simply check the length of requried inputs for the model
    inputs = (key-> km.tensors[key]).(ml.inputs)
    if length(inputs) == 1; out = ml.layer(inputs[1]);
        else; out = ml.layer(inputs...); end
        #else; out = ml.layer(inputs[1]); end

        # SAVE OUTPUTS
        # check if there are multiple outputs (rnn etc.) before saving them to model.tensors
    if length(ml.outputs) == 1; km.tensors[ml.outputs[1]] = out;
        else; for output in ml.outputs; km.tensors[output] = out; end; end
 end

function (m::KnetModel)(args...)

    # resets model.tensors (figure out a smarter/faster reset)
    m.tensors = TensorDict(m.model_layers, m.graph)

    # REGISTER X
    # dumb version

    for (i, arg) in enumerate(args)
        m.tensors[m.model_inputs[i]] = arg
    end

    # LOOP UNTIL ALL TENSORS ARE CALCULATED
    # do until all model.tensors are filled
    # iterate over all layers and call forward on that layer
    while Nothing in values(m.tensors)
        for layer in m.model_layers
            forward(m, layer)
        end
    end

    # Gather outputs
    outs = []
    #m.tensors[m.model_outputs...]
    # DUMB VERSION
    # could be multiple
    if length(m.model_outputs) == 1; outs = m.tensors[m.model_outputs[1]];
        else; for out in m.model_outputs; push!(outs, m.tensors[out]); end;  end

    # return outputs
    return outs

end

"""
    PrintModelTensors(models::KnetModel)
    Displays your model.tensor, showing the size of the tensors that are calculated.
    Might be useful for debugging your model.
"""
function PrintModelTensors(model::KnetModel)
    tensors = model.tensors
    for k in keys(tensors)
        if tensors[k] == Nothing; println(k, "\t=> ", "Nothing")
        else println(k, "\t=> ", size(tensors[k])); end
    end
end

#import Knet: Data

#(m::KnetModel)(x,y) = Knet.nll(m(x), y)

"""
    Calling KnetModel with a Knet.Data object computes the mean nll for
"""
#(m::KnetModel)(d::Data) = Knet.mean(m(x,y) for (x,y) in d)


text:
module Makie

using AbstractPlotting, ImageCore, LinearAlgebra, Statistics, Base64
using GeometryTypes, Colors, ColorVectorSpace, StaticArrays, FixedPointNumbers
import IntervalSets, FileIO
using IntervalSets: ClosedInterval, (..)
using Primes
using Base.Iterators: repeated, drop
using FreeType, FreeTypeAbstraction, UnicodeFun
using PlotUtils, Showoff
using Base: RefValue
import Base: push!, isopen, show

for name in names(AbstractPlotting)
    @eval import AbstractPlotting: $(name)
    @eval export $(name)
end
# Unexported names
using AbstractPlotting: @info, @log_performance, @warn, NativeFont, Key, @key_str

export (..), GLNormalUVMesh
# conflicting identifiers
using AbstractPlotting: Text, volume, VecTypes
using GeometryTypes: widths
export widths, decompose
export VideoStream, recordframe!, finish, record

module ContoursHygiene
    import Contour
end
using .ContoursHygiene
const Contours = ContoursHygiene.Contour

const has_ffmpeg = Ref(false)

function __init__()
    has_ffmpeg[] = try
        success(`ffmpeg -h`)
    catch
        false
    end
end

function logo()
    FileIO.load(joinpath(@__DIR__, "..", "docs", "src", "assets", "logo.png"))
end

include("makie_recipes.jl")
include("utils.jl")
include("glbackend/glbackend.jl")
using .GLMakie
using .GLMakie: assetpath, loadasset

include("video_io.jl")

end


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# ╔═╡ 32465c0e-fcd4-11ea-1544-df26081c7fa7
md"""
# Functions are objects
"""

# ╔═╡ 65c59b14-fcd7-11ea-2a19-3d084b3bca56
square_root = sqrt

# ╔═╡ 249cf7dc-fcdb-11ea-3630-ed2369d20041
square_root(123)

# ╔═╡ 6b91914a-fcd4-11ea-0d27-c99e7ef99354
function double(x)
	x * 2
end

# ╔═╡ b5614d1a-fcd4-11ea-19b9-45043b16b332
function half(x)
	x / 2
end

# ╔═╡ 0991e74c-fce3-11ea-0616-336e1d5d83e9
things = [double, half]

# ╔═╡ 10fe7950-fce3-11ea-1ace-e1676961935e
rand(things)(123)

# ╔═╡ 2424dd62-fce3-11ea-14a6-81792a7dee89
function applyboth(f, g, x)
	f(g(x))
end

# ╔═╡ 2fb8117e-fce3-11ea-2492-55e4768f6e37
applyboth(double, half, 10)

# ╔═╡ a34557ac-fce3-11ea-1391-3d0cddd4201b
md"""
# _map_ and _filter_
"""

# ╔═╡ 70aa854a-fce5-11ea-3477-6df2b0ca1d22
struct Dog
	name
	age
	photo
end

# ╔═╡ cbfc5ede-fce3-11ea-2044-15b8a07ef5f2
data = [
	Dog("Floep", 13, md"![](https://i.imgur.com/4PHFyIE.jpg)"),
	Dog("Hannes", 5, md"![](https://i.imgur.com/nD5c6yF.jpg)"),
	Dog("Fred", 8, md"![](https://i.imgur.com/aYTy1QN.jpg)"),
	Dog("Lily", 3, md"![](https://i.imgur.com/2monPgX.jpg)"),
	Dog("Robert", 15, md"![](https://i.imgur.com/jEwm3Q0.jpg)"),
	Dog("Kit", 2, md"![](https://i.imgur.com/KKtlIEe.jpg)"),
	Dog("Spot", 10, md"![](https://i.imgur.com/hNadBtk.jpg)"),
	]

# ╔═╡ 74f63e2c-fce9-11ea-2145-dd96e9cda96c
md"👉 Show the **photos** of all dogs older than **7 years**."

# ╔═╡ 9b688c4a-fceb-11ea-10b1-590b77c7bfe3
function isold(dog)
	dog.age > 7
end

# ╔═╡ ef6ebf86-fcea-11ea-1118-4f4b4960692b
filter(isold, data)

# ╔═╡ b7608b28-fceb-11ea-3742-a7828971d170
filter(dog -> dog.age > 7, data)

# ╔═╡ c53212da-fceb-11ea-0eeb-617a18323021
special_dogs = filter(data) do dog
	dog.age > 7
end

# ╔═╡ ea0ca73c-fceb-11ea-348a-5df7974b4aba
map(special_dogs) do dog
	dog.photo
end

# ╔═╡ 0758eff0-fcd4-11ea-3186-e1f76a06b91c
bigbreak = html"""

"""

# ╔═╡ 175cb644-fcd5-11ea-22f2-3f96d6d2e637
bigbreak

# ╔═╡ d5c31a1a-fcd8-11ea-0841-1f4a056c048e
bigbreak

# ╔═╡ c8fc9460-fce3-11ea-0f2a-1b87abdd12b7
bigbreak

# ╔═╡ Cell order:
# ╟─32465c0e-fcd4-11ea-1544-df26081c7fa7
# ╠═65c59b14-fcd7-11ea-2a19-3d084b3bca56
# ╠═249cf7dc-fcdb-11ea-3630-ed2369d20041
# ╟─175cb644-fcd5-11ea-22f2-3f96d6d2e637
# ╠═6b91914a-fcd4-11ea-0d27-c99e7ef99354
# ╠═b5614d1a-fcd4-11ea-19b9-45043b16b332
# ╠═0991e74c-fce3-11ea-0616-336e1d5d83e9
# ╠═10fe7950-fce3-11ea-1ace-e1676961935e
# ╠═2424dd62-fce3-11ea-14a6-81792a7dee89
# ╠═2fb8117e-fce3-11ea-2492-55e4768f6e37
# ╟─d5c31a1a-fcd8-11ea-0841-1f4a056c048e
# ╟─a34557ac-fce3-11ea-1391-3d0cddd4201b
# ╠═70aa854a-fce5-11ea-3477-6df2b0ca1d22
# ╟─cbfc5ede-fce3-11ea-2044-15b8a07ef5f2
# ╟─74f63e2c-fce9-11ea-2145-dd96e9cda96c
# ╠═9b688c4a-fceb-11ea-10b1-590b77c7bfe3
# ╠═ef6ebf86-fcea-11ea-1118-4f4b4960692b
# ╠═b7608b28-fceb-11ea-3742-a7828971d170
# ╠═c53212da-fceb-11ea-0eeb-617a18323021
# ╠═ea0ca73c-fceb-11ea-348a-5df7974b4aba
# ╟─c8fc9460-fce3-11ea-0f2a-1b87abdd12b7
# ╟─0758eff0-fcd4-11ea-3186-e1f76a06b91c

max_stars_repo_name: mathijsvdv/ComputationalThinking

max_issues_repo_name: mathijsvdv/ComputationalThinking

max_forks_repo_name: mathijsvdv/ComputationalThinking


text:
using ArgParse
s = ArgParseSettings()
@add_arg_table s begin
    "--part2"
    help = "Do part 2"
    action = :store_true
    "--test"
    help = "Run test data"
    action = :store_true
end
parsed_args = parse_args(ARGS, s)

test_suffix = parsed_args["test"] ? "_test" : ""

points = Set()
axes = []

for l in readlines("input/day_13$test_suffix.txt")
    if contains(l, ",")
        p = split(l, ",")
        push!(points, (x = parse(UInt16, p[1]), y = parse(UInt16, p[2])))
    elseif length(l) > 4
        axis = l[begin+11:end]
        p = split(axis, "=")
        push!(axes, (axis = p[1], line = parse(UInt16, p[2])))
    end
end

function dofold(points, fold)
    result = Set()
    for p in points
        if fold.axis == "x"
            x = p.x
            if x > fold.line
                x = 2 * fold.line - x
            end
            push!(result, (x = x, y = p.y))
        else
            y = p.y
            if y > fold.line
                y = 2 * fold.line - y
            end
            push!(result, (x = p.x, y = y))
        end
    end

    result
end

function part1()
    length(dofold(points, axes[1]))
end

function part2()
    p = points
    for fold in axes
        p = dofold(p, fold)
    end
    mx = 0
    my = 0
    for pt in p
        mx = max(pt.x, mx)
        my = max(pt.y, my)
    end
    code = zeros(mx + 1, my + 1)
    for pt in p
        code[pt.x+1, pt.y+1] = 1
    end
    for y = 1:my+1
        line = Vector{Char}()
        for x = 1:mx+1
            push!(line, code[x, y] == 1 ? '@' : ' ')
        end
        println(String(line))
    end
end

p = parsed_args["part2"] ? part2 : part1
answer = p()
println("The answer is: $answer")


text:
using Test
import Literate

const EXAMPLES_DIR = @__DIR__
const OUTPUT_DIR   = joinpath(@__DIR__, "generated")

const EXAMPLES = [
    "Convex hull and intersection.jl",
    "Extended Formulation.jl",
    "Minimal Robust Positively Invariant Set.jl"
]

@testset "test_examples.jl" begin
    @testset "$(example)" for example in EXAMPLES
        Literate.script(joinpath(EXAMPLES_DIR, example), OUTPUT_DIR)
        include(joinpath(OUTPUT_DIR, example))
    end
end

max_stars_repo_name: blegat/Polyhedra

max_issues_repo_name: blegat/Polyhedra

max_forks_repo_name: blegat/Polyhedra


text:
# push!(LOAD_PATH,"../src/")

using Documenter, OpenQuantumSystems

makedocs(
    modules = [OpenQuantumSystems],
    format = Documenter.HTML(
        prettyurls = !("local" in ARGS),
        canonical = "https://detrin.github.io/OpenQuantumSystems.jl/latest/",
        analytics = "G-F12F2Y92VY",
    ),
    clean = false,
    sitename = "OpenQuantumSystems.jl",
    authors = "Daniel Herman",
    linkcheck = !("skiplinks" in ARGS),
    pages = [
        "Home" => "index.md",
        "Documentation" => "documentation.md",
        "Tutorials" => [
            "Dimer" => "tutorials/dimer.md"
            ]
    ],
)

deploydocs(repo = "github.com/detrin/OpenQuantumSystems.jl.git")


text:
module TestExamples3DStructured

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "structured_3d_dgsem")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

@testset "Structured mesh" begin
  @trixi_testset "elixir_advection_basic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"),
      # Expected errors are exactly the same as with TreeMesh!
      l2   = [0.00016263963870641478],
      linf = [0.0014537194925779984])
  end

  @trixi_testset "elixir_advection_free_stream.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_free_stream.jl"),
      l2   = [1.2908196366970896e-14],
      linf = [1.0262901639634947e-12],
      atol = 8e-13, # required to make tests pass on Windows
      )
  end

  @trixi_testset "elixir_advection_nonperiodic_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonperiodic_curved.jl"),
      l2   = [0.0004483892474201268],
      linf = [0.009201820593762955])
  end

  @trixi_testset "elixir_advection_restart.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"),
      l2   = [0.0025903889347585777],
      linf = [0.018407576968841655])
  end

  @trixi_testset "elixir_euler_source_terms.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
      # Expected errors are exactly the same as with TreeMesh!
      l2   = [0.010385936842224346, 0.009776048833895767, 0.00977604883389591, 0.009776048833895733, 0.01506687097416608],
      linf = [0.03285848350791731, 0.0321792316408982, 0.032179231640894645, 0.032179231640895534, 0.0655408023333299])
  end

  @trixi_testset "elixir_euler_free_stream.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream.jl"),
      l2   = [2.8815700334367128e-15, 9.361915278236651e-15, 9.95614203619935e-15, 1.6809941842374106e-14, 1.4815037041566735e-14],
      linf = [4.1300296516055823e-14, 2.0444756998472258e-13, 1.0133560657266116e-13, 2.0627943797535409e-13, 2.8954616482224083e-13])
  end

  @trixi_testset "elixir_euler_free_stream.jl with FluxRotated(flux_lax_friedrichs)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream.jl"),
      surface_flux=FluxRotated(flux_lax_friedrichs),
      l2   = [2.8815700334367128e-15, 9.361915278236651e-15, 9.95614203619935e-15, 1.6809941842374106e-14, 1.4815037041566735e-14],
      linf = [4.1300296516055823e-14, 2.0444756998472258e-13, 1.0133560657266116e-13, 2.0627943797535409e-13, 2.8954616482224083e-13])
  end

  @trixi_testset "elixir_euler_source_terms_nonperiodic_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonperiodic_curved.jl"),
    l2   = [0.0032940531178824463, 0.003275679548217804, 0.0030020672748714084, 0.00324007343451744, 0.005721986362580164],
    linf = [0.03156756290660656, 0.033597629023726316, 0.02095783702361409, 0.03353574465232212, 0.05873635745032857])
  end

  @trixi_testset "elixir_euler_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.011367083018614027, 0.007022020327490176, 0.006759580335962235, 0.006820337637760632, 0.02912659127566544],
      linf = [0.2761764220925329, 0.20286331858055706, 0.18763944865434593, 0.19313636558790004, 0.707563913727584],
      tspan = (0.0, 0.25),
      coverage_override = (polydeg=3,)) # Prevent long compile time in CI
  end

  @trixi_testset "elixir_euler_sedov.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"),
      l2   = [5.30310390e-02, 2.53167260e-02, 2.64276438e-02, 2.52195992e-02, 3.56830295e-01],
      linf = [6.16356950e-01, 2.50600049e-01, 2.74796377e-01, 2.46448217e-01, 4.77888479e+00],
      tspan = (0.0, 0.3))
  end

  @trixi_testset "elixir_mhd_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"),
      l2   = [0.009082353036644902, 0.007128360240528109, 0.006970330025996491, 0.006898850266874514,
              0.03302008823756457, 0.003203389099143526, 0.003077498677885352, 0.0030740006760477624,
              4.192129696970217e-5],
      linf = [0.2883946030582689, 0.25956437344015054, 0.2614364943543665, 0.24617277938134657,
              1.1370443512475847, 0.1278041831463388, 0.13347391885068594, 0.1457563463643099,
              0.0021174246048172563],
      tspan = (0.0, 0.25))
  end

  @trixi_testset "elixir_mhd_alfven_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [0.003015476175153681, 0.00145499403283373, 0.0009125744757935803, 0.0017703080480578979,
              0.0013046447673965966, 0.0014564863387645508, 0.0013332311430907598, 0.001647832598455728,
              0.0013647609788548722],
      linf = [0.027510637768610846, 0.02797062834945721, 0.01274249949295704, 0.038940694415543736,
              0.02200825678588325, 0.03167600959583505, 0.021420957993862344, 0.03386589835999665,
              0.01888303191983353],
      # Use same polydeg as everything else to prevent long compile times in CI
      coverage_override = (polydeg=3,))
  end

  @trixi_testset "elixir_mhd_alfven_wave.jl with flux_lax_friedrichs" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [0.003047854479955232, 0.0014572199588782184, 0.0009093737183251411, 0.0017937548694553895,
              0.0013010437110755424, 0.0014545607744895874, 0.001328514015121245, 0.001671342529206066,
              0.0013653963058149186],
      linf = [0.027719103797310463, 0.027570111789910784, 0.012561901006903103, 0.03903568568480584,
              0.021311996934554767, 0.03154849824135775, 0.020996033645485412, 0.03403185137382961,
              0.019488952445771597],
      surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell),
      # Use same polydeg as everything else to prevent long compile times in CI
      coverage_override = (polydeg=3,))
  end

  @trixi_testset "elixir_mhd_ec_shockcapturing.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec_shockcapturing.jl"),
      l2   = [0.009352631220872144, 0.008058649103542618, 0.008027041293333663, 0.008071417851552725,
              0.034909149665869485, 0.00393019428600812, 0.0039219074393817, 0.003906321245184237,
              4.197255300781248e-5],
      linf = [0.30749098250807516, 0.2679008863509767, 0.271243087484388, 0.26545396569129537,
              0.9620950892188596, 0.18163281157498123, 0.15995708312378454, 0.17918221526906408,
              0.015138346608166353],
      tspan = (0.0, 0.25),
      # Use same polydeg as everything else to prevent long compile times in CI
      coverage_override = (polydeg=3,))
  end
end

# Clean up afterwards: delete Trixi output directory
@test_nowarn rm(outdir, recursive=true)

end # module


text:
include("LinearAlgebra/FakeFmpqMat.jl")
include("LinearAlgebra/LatEnum.jl")
include("LinearAlgebra/Howell.jl")
include("LinearAlgebra/Sparse.jl")
#include("LinearAlgebra/AbGrp.jl")


text:
##############################################################################
#
# REFERENCES: "Statistical Distributions"
#
##############################################################################

struct EmpiricalUnivariateDistribution <: ContinuousUnivariateDistribution
	values::Vector{Float64}
    support::Vector{Float64}
	cdf::Function
	entropy::Float64
	kurtosis::Float64
	mean::Float64
	median::Float64
	modes::Vector{Float64}
	skewness::Float64
	var::Float64
end

@distr_support EmpiricalUnivariateDistribution d.values[1] d.values[end]

function EmpiricalUnivariateDistribution(x::Vector)
    sx = sort(x)
	EmpiricalUnivariateDistribution(sx,
                                    unique(sx),
	                                ecdf(x),
	                                NaN,
	                                NaN,
	                                mean(x),
	                                median(x),
	                                Float64[],
	                                NaN,
	                                var(x))
end

entropy(d::EmpiricalUnivariateDistribution) = d.entropy

kurtosis(d::EmpiricalUnivariateDistribution) = d.kurtosis

mean(d::EmpiricalUnivariateDistribution) = d.mean

median(d::EmpiricalUnivariateDistribution) = d.median

modes(d::EmpiricalUnivariateDistribution) = Float64[]

skewness(d::EmpiricalUnivariateDistribution) = NaN

var(d::EmpiricalUnivariateDistribution) = d.var

### Evaluation

cdf(d::EmpiricalUnivariateDistribution, x::Float64) = d.cdf(x)

function pdf(d::EmpiricalUnivariateDistribution, x::Float64)
    ## TODO: Create lookup table for discrete case
    1.0 / length(d.values)
end

function quantile(d::EmpiricalUnivariateDistribution, p::Float64)
    n = length(d.values)
    index = floor(Int,p * n) + 1
    index > n ? d.values[n] : d.values[index]
end

function rand(d::EmpiricalUnivariateDistribution)
    d.values[rand(1:length(d.values))]
end

### fit model

function fit_mle(::Type{EmpiricalUnivariateDistribution},
             x::Vector{T}) where T <: Real
	EmpiricalUnivariateDistribution(x)
end


text:
module MPStates

using LinearAlgebra, TensorOperations, Random, HDF5, Printf, KrylovKit

include("./cache.jl")
include("./tensor_contractions.jl")
include("./tensor_factorizations.jl")
include("./mps.jl")
include("./mps_operations.jl")
include("./mpo.jl")
include("./mpo_operations.jl")
include("./dmrg_opts.jl")
include("./dmrg.jl")
include("./hamiltonians.jl")

end # module MPStates


text:
## test.jl
## (c) 2013--2016 David A. van Leeuwen

## various tests for NamedArrays

## This code is licensed under the MIT license
## See the file LICENSE.md in this distribution

## NamedArrays is loaded by runtests.jl

print("Starting test, no assertions should fail... ")

include("bugs.jl")

include("base.jl")

include("constructors.jl")

include("arithmetic.jl")

include("convert.jl")

include("index.jl")

include("names.jl")

include("keepnames.jl")

print("copy, ")
## copy
m = copy(n)
@test m == n
copyto!(m, n)
@test m == n
m = deepcopy(n)
@test m == n

print("sum, ")
## sum
@test sum(n) == sum(n.array)
@test sum(n, dims=1).array == sum(n.array, dims=1)
@test sum(n, dims=2).array == sum(n.array, dims=2)
@test names(sum(n, dims=2), 1) == ["one", "two"]

println("conversions")
## conversion
@test convert(Array, n) == n.array
@test map(Float32, n).array == map(Float32, n.array)

println("changing names: ")
## changingnames
for  f = (:sum, :prod, :maximum, :minimum, :mean, :std, :var)
    println("    ", f)
    for dim in 1:2
        @eval @test ($f)(n.array, dims=$dim) == ($f)(n, dims=$dim).array
        @eval @inferred ($f)(n, dims=$dim)
    end
end

for f in (:cumprod, :cumsum)
    for dim in 1:2
        @eval @test ($f)(n.array, dims=$dim) == ($f)(n, dims=$dim).array
        @eval @inferred ($f)(n, dims=$dim)
    end
end

# NOTE: KahanSummation do not support Julia 0.7 dims keyword argument at the moment
@test cumsum_kbn(n.array, 1) == cumsum_kbn(n, 1).array
@test cumsum_kbn(n.array, 2) == cumsum_kbn(n, 2).array
@inferred cumsum_kbn(n, 1)
@inferred cumsum_kbn(n, 2)

print("multi-dimensional, ")
#multidimensional
m = NamedArray(rand(2,3,4,3,2,3))
for i1 in 1:2, i2 in 1:3, i3 in 1:4, i4 in 1:3, i5 in 1:2, i6 in 1:3
    @test m[i1,i2,i3,i4,i5,i6] == m.array[i1,i2,i3,i4,i5,i6]
    @test m[string(i1), string(i2), string(i3), string(i4), string(i5), string(i6)] == m.array[i1,i2,i3,i4,i5,i6]
end
m[1, :, 2, :, 2, 3].array == m.array[1, :, 2, :, 2, 3]
i = [3 2 4; 1 4 2]
@test n[:, i].array == n.array[:, i]
@test names(n[:, i], 1) == names(n, 1)
@test names(n[:, i], 2) == map(string, 1:2)
@test names(n[:, i], 3) == map(string, 1:3)

print("dodgy indices, ")
## weird indices
m = NamedArray(rand(4), ([1//1, 1//2, 1//3, 1//4],), ("rational",))
@test m[1//2] == m.array[2]
@test m[[1//4,1//3]] == m.array[[4,3]]
m[1//4] = 1
@test m[4] == 1

m = NamedArray(rand(4), ([4, 3, 2, 1],), ("reverse confusion",))
@test m[1] == m.array[1] ## Integer precedence
@test m[Name(1)] == m.array[4] ## but Name() indicates names
## this goes wrong for julia-v0.3
@test m[[4,3,2,1]].array == m.array[[4, 3, 2, 1]]
@test m[Name.([4, 3, 2, 1])].array == m

print("sort, ")
m = NamedArray(rand(100))
for rev in [false, true]
    ms = sort(m, rev=rev)
    @test ms.array == sort(m.array, rev=rev)
    @test names(ms, 1) == names(m, 1)[sortperm(m.array, rev=rev)]
end
m = NamedArray(rand(10,10))
for rev in [false, true]
    for dim in 1:2
        ms = sort(m, dims=dim, rev=rev)
        @test ms.array == sort(m.array, dims=dim, rev=rev)
        @test names(ms, dim) == [string(i) for i in 1:size(ms, dim)]
    end
end

print("broadcast, ")
@test broadcast(-, n, mean(n, dims=1)).array == broadcast(-, n.array, mean(n.array, dims=1))

print("vectorized, ")

## a selection of vectorized functions
for f in  (:sin, :cos, :tan,  :sinpi, :cospi, :sinh, :cosh, :tanh, :asin, :acos, :atan,
           :sinc, :cosc, :deg2rad, :log, :log2, :log10, :log1p, :exp, :exp2, :exp10,
           :expm1, :abs, :abs2, :sign, :sqrt)
    @eval begin
        m = ($f).(n)
        @test m.array == ($f).(n.array)
        @test namesanddim(m) == namesanddim(n)
    end
end

include("rearrange.jl")

print("eachindex, ")
## eachindex
for I in eachindex(n)
    @test n[I] == n.array[I]
end

include("matrixops.jl")

include("show.jl")

include("speed.jl")

# julia issue #17328
a = NamedArray([1.0, 2.0, 3.0, 4.0])
@test sum(abs, a, dims=1)[1] == 10

println("done!")


text:
# ---
# title: Custom Title
# ---

# This is the content


text:
"""
    DistanceGradient

The `sources` field is the *linear* indices of the matrix, from which the
distance must be calculated.
"""
struct DistanceGradient <: NeutralLandscapeMaker
    sources::Vector{Integer}
end

DistanceGradient() = DistanceGradient([1])

function _landscape!(mat, alg::DistanceGradient)
    @assert maximum(alg.sources) <= length(mat)
    @assert minimum(alg.sources) > 0
    coordinates = Matrix{Float64}(undef, (2, prod(size(mat))))
    for (i, p) in enumerate(Iterators.product(axes(mat)...))
        coordinates[1:2, i] .= p
    end
    tree = KDTree(coordinates[:,alg.sources])
    mat[:] .= nn(tree, coordinates)[2]
    return mat
end


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

module PkgTestSuite

using Coverage
using Documenter
using TimerOutputs
if VERSION < v"1.0.0"
    using Base.Test
else
    using Test
end

USING_LINT = true
try
    using Lint
catch
    info("Cannot use Lint.jl with Julia version $VERSION")
    USING_LINT = false
end

include("checkheader.jl")
include("checktabs.jl")

"""
    determine_pkg_name(pkg::String="")

A helper function to determine package name.
"""
function determine_pkg_name(pkg::String="")
    if pkg == ""
        if !haskey(ENV, "TRAVIS_REPO_SLUG")
            error("package name not given and not in Travis-CI environmnt")
        end
        s = ENV["TRAVIS_REPO_SLUG"]
        info("TRAVIS_REPO_SLUG = $s")
        pkg = replace(split(s,'/')[end], ".jl", "")
    end
    return pkg
end

"""
    init()

Clone and build package subject to testing.
"""
function init(pkg_name::String="")
    pkg_name = determine_pkg_name(pkg_name)
    info("PkgTestSuite.init(): determined package name to be $pkg_name")
    inside_travis = haskey(ENV, "TRAVIS")
    if inside_travis
        pkg_path = ENV["TRAVIS_BUILD_DIR"]
    else
        pkg_path = Pkg.dir(pkg_name)
    end
    info("PkgTestSuite.init(): detemined location of package to be $pkg_path")
    if pkg_name == "PkgTestSuite"
        # self-checking PkgTestSuite, not build another time even inside Travis-CI
        return nothing
    end
    repo = LibGit2.GitRepo(Pkg.dir("PkgTestSuite"))
    hash_long = string(LibGit2.revparseid(repo, "master"))
    hash_short = hash_long[1:7]
    info("PkgTestSuite.init(): Determined PkgTestSuite commit hash to be $hash_short")
    if inside_travis
        info("Inside Travis-CI, cloning and building package")
        Pkg.clone(pkg_path)
        Pkg.build(pkg_name)
    end
    return nothing
end

"""
    test(pkg_name)

Run tests for package.
"""
function test(pkg_name=""; run_tests=true)

    pkg_name = determine_pkg_name(pkg_name)
    pkg_dir = Pkg.dir(pkg_name)
    cd(pkg_dir)
    if VERSION < v"0.7.0-beta2.0"
        Base.require(Symbol(pkg_name))
        pkg = getfield(Main, Symbol(pkg_name))
    else
        pkg = Base.require(Module(Symbol(pkg_name)), Symbol(pkg_name))
    end

    # check for headers and tabulators
    checkheader(pkg_name)
    checktabs(pkg_name)

    # make it possible to allow lint check or doctest fail
    # without failing whole build
    strict_lint = (get(ENV, "LINT_STRICT", "true") == "true")
    strict_docs = (get(ENV, "DOCUMENTER_STRICT", "true") == "true")

    # run lint
    if USING_LINT
        results = lintpkg(pkg_name)
        if !isempty(results)
            info("Lint.jl is a tool that uses static analysis to assist ",
                 "in the development process by detecting common bugs and ",
                 "potential issues.")
            info("For this package, Lint.jl report is following:")
            display(results)
            info("For more information, see https://lintjl.readthedocs.io/en/stable/")
            warn("Package syntax test has failed.")
            if strict_lint
                errors_and_warnings = filter(i -> !isinfo(i), results)
                @test isempty(errors_and_warnings)
            end
        else
            info("Lint.jl: syntax check pass.")
        end
    else
        info("Lint.jl is not $VERSION compatible yet.")
    end

    # generate documentation and run doctests
    docs_dir = Pkg.dir(pkg_name, "docs")
    docs_src_dir = Pkg.dir(pkg_name, "docs", "src")
    readme_file = Pkg.dir(pkg_name, "README.md")
    index_file = Pkg.dir(pkg_name, "docs", "src", "index.md")
    if !isdir(docs_dir)
        info("Creating new directory $docs_dir")
        mkdir(docs_dir)
    end
    if !isdir(docs_src_dir)
        info("Creating new directory $docs_src_dir")
        mkdir(docs_src_dir)
    end
    if !isfile(index_file)
        warn("$index_file not found.")
        if isfile(readme_file)
            info("Copying $readme_file to $index_file")
            cp(readme_file, index_file)
        else
            index_file_default_contents = """# $pkg_name.jl

Package documentation missing. Start writing documentation to package by
creating docs/src/index.md where it is described what this package does.

Also, create README.md"""
            open(index_file, "w") do fid
                write(fid, index_file_default_contents)
            end
        end
    end
    cd(docs_dir)
    makefile = joinpath(docs_dir, "make.jl")
    if isfile(makefile)
        include(makefile)
    else
        warn("$makefile not found, using defaults to generate documentation.")
        makedocs(
            modules = [pkg],
            format = :html,
            checkdocs = :all,
            sitename = "$pkg_name.jl",
            pages = ["index.md"],
            strict = strict_docs)
    end

    # run pkg tests
    runtests_file = Pkg.dir(pkg_name, "test", "runtests.jl")
    if isfile(runtests_file) && run_tests
        Pkg.test(pkg_name, coverage=true)
    else
        warn("Not running tests for $pkg_name.")
    end
    return nothing
end

"""
    deploy(pkg::String)

Deploy package.
"""
function deploy(pkg::String="")

    pkg = determine_pkg_name(pkg)

    if !haskey(ENV, "TRAVIS")
        info("Looks that you are not running deploy on CI platform.")
        info("Documentation generated to ", Pkg.dir(pkg, "docs", "site"))
        return nothing
    end

    # upload results to coveralls.io
    cd(Pkg.dir(pkg))
    result = Coveralls.process_folder()
    Coveralls.submit(result)

    # deploy documentation to juliafem.github.io
    cd(Pkg.dir(pkg, "docs"))
    deploy_file = Pkg.dir(pkg, "docs", "deploy.jl")
    if isfile(deploy_file)
        include(deploy_file)
    else
        warn("deploy.jl not found, deploying using default settings")
        deploydocs(
            repo = "github.com/JuliaFEM/$pkg.jl.git",
            julia = "0.6",
            deps = nothing,
            make = nothing)
    end

    return nothing
end

export init, test, deploy

end


text:
#gnum-constants.jl
#setting gnum values to unum constants.

function nan!{ESS,FSS}(x::Gnum{ESS,FSS})
  x.scratchpad.flags |= GNUM_NAN_MASK
end

@generated function inf!{ESS,FSS,side}(x::Gnum{ESS,FSS}, flags::UInt16, ::Type{Val{side}})
  :(clear_gflags!(x.$side); x.$side.flags = flags | GNUM_INF_MASK; nothing)
end
@generated function mmr!{ESS,FSS,side}(x::Gnum{ESS,FSS}, flags::UInt16, ::Type{Val{side}})
  :(clear_gflags!(x.$side); x.$side.flags = flags | GNUM_MMR_MASK; nothing)
end
@generated function sss!{ESS,FSS,side}(x::Gnum{ESS,FSS}, flags::UInt16, ::Type{Val{side}})
  :(clear_gflags!(x.$side); x.$side.flags = flags | GNUM_SSS_MASK; nothing)
end
@generated function zero!{ESS,FSS,side}(x::Gnum{ESS,FSS}, ::Type{Val{side}})
  :(clear_gflags!(x.$side); x.$side.flags = GNUM_ZERO_MASK; nothing)
end

#testing gnum values for unum constants
function is_nan{ESS,FSS}(x::Gnum{ESS,FSS})
  x.scratchpad.flags & GNUM_NAN_MASK != 0
end

set_g_mmr!{ESS,FSS}(x::Unum{ESS,FSS}) = (x.flags |= GNUM_MMR_MASK)

is_g_inf{ESS,FSS}(x::Unum{ESS,FSS}) = x.flags & GNUM_INF_MASK != 0
@generated function is_inf{ESS,FSS,side}(x::Gnum{ESS,FSS}, ::Type{Val{side}})
 :(is_g_inf(x.$side))
end

is_g_mmr{ESS,FSS}(x::Unum{ESS,FSS}) = x.flags & GNUM_MMR_MASK != 0
@generated function is_mmr{ESS,FSS,side}(x::Gnum{ESS,FSS}, ::Type{Val{side}})
  :(is_g_mmr(x.$side))
end

is_g_sss{ESS,FSS}(x::Unum{ESS,FSS}) = x.flags & GNUM_SSS_MASK != 0
@generated function is_sss{ESS,FSS,side}(x::Gnum{ESS,FSS}, ::Type{Val{side}})
  :(is_g_sss(x.$side))
end

is_g_zero{ESS,FSS}(x::Unum{ESS,FSS}) = x.flags & GNUM_ZERO_MASK != 0
@generated function is_zero{ESS,FSS,side}(x::Gnum{ESS,FSS}, ::Type{Val{side}})
  :(is_g_zero(x.$side))
end

@generated function force_from_flags!{ESS,FSS,side}(src::Gnum{ESS,FSS}, dest::Unum{ESS,FSS}, ::Type{Val{side}})
  quote
    (src.$side.flags & GNUM_INF_MASK != 0) && (inf!(dest, src.$side.flags & UNUM_SIGN_MASK); return true)
    (src.$side.flags & GNUM_MMR_MASK != 0) && (mmr!(dest, src.$side.flags & UNUM_SIGN_MASK); return true)
    (src.$side.flags & GNUM_SSS_MASK != 0) && (sss!(dest, src.$side.flags & UNUM_SIGN_MASK); return true)
    (src.$side.flags & GNUM_ZERO_MASK != 0) && (zero!(dest); return true)
    false
  end
end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl