text:
function ctns2cmps(V,α,χ)
    Q = zeros(Float64, χ, χ)
    R = zeros(Float64, χ, χ)
    s = (V^2 - abs(α)^4)^(1/4)
    for n = 1:χ
        Q[n,n] += 4*(n-1)^2*sqrt(s)
        Q[n,n] -= sqrt(2)*2*(n-1)*sqrt(n-1)
        Q[n,n] += 1/(2*s)*(V+1)*(2*n-1)
        if n+1 <= χ
            R[n+1,n] = α/s*sqrt(n/2)
        end
        if n-1 > 0
            R[n-1,n] = α/s*sqrt((n-1)/2)
        end
        if n+2 <= χ
            Q[n+2,n] -= sqrt(2)*(n-1)*sqrt(n)
            Q[n+2,n] += 1/(2*s)*(V+1)*sqrt((n+1)*(n))
        end
        if n-2 > 0
            Q[n-2,n] -= sqrt(2)*(n-3)*sqrt(n-2)
            Q[n-2,n] += 1/(2*s)*(V+1)*sqrt((n-1)*(n-2))
        end
    end
    return Q, R
end

function ctns2cmps(V,α,χ)
    M = real([V (- α * adjoint(α) / 2); (- conj(α) * transpose(α) / 2) conj(V)])
    Msq = sqrt(M)
    if !isreal(Msq)
        error("invalid state")
    end
    Ω = real(Msq)

    D,U = eigen(M)

    Q = zeros(Float64, χ, χ)
    R = zeros(Float64, χ, χ)
    s = (V^2 - abs(α)^4)^(1/4)
    for n = 1:χ
        Q[n,n] += 4*(n-1)^2*sqrt(s)
        Q[n,n] -= sqrt(2)*2*(n-1)*sqrt(n-1)
        Q[n,n] += 1/(2*s)*(V+1)*(2*n-1)
        if n+1 <= χ
            R[n+1,n] = α/s*sqrt(n/2)
        end
        if n-1 > 0
            R[n-1,n] = α/s*sqrt((n-1)/2)
        end
        if n+2 <= χ
            Q[n+2,n] -= sqrt(2)*(n-1)*sqrt(n)
            Q[n+2,n] += 1/(2*s)*(V+1)*sqrt((n+1)*(n))
        end
        if n-2 > 0
            Q[n-2,n] -= sqrt(2)*(n-3)*sqrt(n-2)
            Q[n-2,n] += 1/(2*s)*(V+1)*sqrt((n-1)*(n-2))
        end
    end
    return Q, R
end


text:
using SpecialMatrices
using Test
using LinearAlgebra

const files = (
    "companion",
    "frobenius",
    "hankel",
    "hilbert",
    "strang",
    "toeplitz",
    "vandermonde",
)

@testset "SpecialMatrices.jl" begin
    @testset "$(titlecase(f)) matrix" for f in files
        include("$f.jl")
    end
end

ext: jl


text:
@testset "MochaRegression" begin
    using Mocha
    using SGMCMC
    using MLUtilities
    using MochaRegressionWrapper
    using MochaRegression
    using MochaModelFactories
    n = 100 # size of training sample
    nn = 100
    x = randn(10,n+nn)
    y = reshape(sin(x[1,:]+x[2,:])+cos(x[3,:]+x[4,:]),(1,n+nn))  # two outputs, only first 4 inputs are relevant
    y = y + 0.5*randn(size(y))                      # add noise: a good fit would give RMSE=0.5
    # split the training and testing samples
    ytest = y[:,1:nn]      # test
    xtest = x[:,1:nn]      # test
    ytrain = y[:,nn+1:end]   # training
    xtrain= x[:,nn+1:end]    # training
    backend = MochaRegressionWrapper.initMochaBackend(false)
    nunits = 10
    model,name = make_dense_nn_regression([nunits],1)
    dm = MochaRegressionModel(xtrain,ytrain,model,backend)
    dmtest = MochaRegressionModel(xtest,ytest,model,backend,do_shuffle=false)#  test set accuracy
    dmteststorage = MochaRegressionModel(xtest,ytest,model,backend,do_shuffle=false,do_storage=true)
    # initialised from models
    @test 10*nunits+nunits+nunits*1+1 == MochaRegression.fetchnparams(dm)

    # test setparams! and getparams
    d = MochaRegression.fetchnparams(dm)
    x = randn(d)
    MochaRegression.setparams!(dm,x)
    @test MochaRegression.getparams(dm) == x
    g = DataModel.getgrad(dm)
    l = DataModel.getllik(dm)
    for i in 1:10
        x = 0.01*randn(d)
        @test abs(checkgrad(x,l,g)) < 1e-2
    end

    # test evaluate functions
    predictions = MochaRegression.getpredictions(dmteststorage,x)
    manual_results = MochaRegression.evaluate(dmteststorage,predictions[:predictions])
    mocha_results = MochaRegression.evaluate(dmtest,x)
    @test abs(manual_results[:loss] -mocha_results[:loss]) < 1
end


text:
module FIEREE

using PyCall, JSON, Statistics, LinearAlgebra, DataFrames, Dates, HTTP, ProgressMeter, Polynomials

ee = pyimport("ee")

struct Domain
    bbox::Vector
    resolution::Float64
    shape::Vector
    crs::String
end

function round_out(bbox::Vector,resolution::Float64)
    min_buff = bbox .- (bbox .% resolution)
    min_buff = map(x->round.(x,digits=6), min_buff)

    max_buff = bbox .+ (resolution .- (bbox .% resolution))
    max_buff = map(x->round.(x,digits=6), max_buff)

    new_bbox = vcat(min_buff[1:2],max_buff[3:end])
end

function get_domain(bbox::Vector,resolution::Float64; crs::String="EPSG:4326")
    if any((bbox .% resolution) .> 0 )
        bbox = round_out(bbox,resolution)
    end

    xdim = convert(Int,round((bbox[3]-bbox[1]) / resolution))
    ydim = convert(Int,round((bbox[4]-bbox[2]) / resolution))
    shape = [xdim,ydim]

    Domain(bbox,resolution,shape,crs)
end

function domain_to_coords(domain::Domain)
    res = domain.resolution
    west,east = domain.bbox[1], (domain.bbox[3] - res)
    south,north = (domain.bbox[2] - res), domain.bbox[4]
    xx = collect(Array{Float32}(west:res:east))
    yy = collect(Array{Float32}(south:res:north))
    return xx,yy
end

function get_ee_region(bbox::Vector)
    minx, miny, maxx, maxy = bbox
    coordinates = [minx,miny,minx,maxy,maxx,maxy,maxx,miny,minx,miny]

    ee.Geometry.Polygon(coordinates...)
end

function geom_to_bbox(geom)
    coords = geom.geometry().bounds().coordinates().getInfo()[1,:,:]
    minx, maxx = min(coords[:,1]...), max(coords[:,1]...)
    miny, maxy = min(coords[:,2]...), max(coords[:,2]...)

    return [minx,miny,maxx,maxy]

end

function response_to_array(response,width,height;channels=1,dtype=Float64)
    buffer = IOBuffer()
    write(buffer,response.content)
    foo = reinterpret(dtype,take!(buffer))
    close(buffer) # close the io buffer for clean memory management

    # find the number of pixels we expect based on img
    data_size::Int64 = width*height*channels
    # find the offset of what was returned vs what is expected
    offset::Int64 = size(foo,1)-data_size + 1

    arr::Array{dtype} = dtype.(foo)[offset:end]

    if channels > 1
        # reshape array from 1-d to correct 3-d dimensions
        arr3d::Array{dtype} = reshape(arr,(channels,height,width))
        # reorder the dimensions for row-major to column-major
        img = permutedims(arr3d,(1,3,2))
    else
        arr2d::Array{dtype} = reshape(arr,(width,height))
        img = rotl90(arr2d)
    end

    # cast the result as a float32 array
    img = convert(Array{Float32}, img)

end

function exponential_backoff(c::Int)
    N = 2^c -1
    backoffTime = (1/(N+1)) * sum(Array(1:1:N))
    bt = max(backoffTime,0.1) # force backoffs to be at least 0.1s
end

function get_s1_dates(domain::Domain,start_time::String,end_time::String)
    geom = get_ee_region(domain.bbox)
    s1 = (ee.ImageCollection("COPERNICUS/S1_GRD").
        filterDate(start_time,end_time).
        filterBounds(geom).
        filter(ee.Filter.eq("orbitProperties_pass","ASCENDING")).
        select("VV").
        sort("system:time_start")
    )

    date_str = ee.List(s1.map(img_to_date).aggregate_array("isodate")).distinct().getInfo()
    # convert epoch milliseconds to date
    # Julia considers start of epoch Jan 1, 0000 so we add 1970 years to convert to correct times
    all_dates = map(x->DateTime(x),date_str)

    expected_dates = collect(all_dates[1]:Dates.Day(12):all_dates[end])

    mask = sum(map(x -> x.==expected_dates,all_dates),dims=1)[1]
    indices = findall(!iszero,mask)

    final_dates = expected_dates[indices]
end

# helper function to set date info, used for getting list of dates
function img_to_date(img)
    return img.set("isodate",img.date().format("YYYY-MM-dd"))
end

# helper function to mask high view angle areas
function s1_qa(img)
    angle = img.select("angle")
    return img.updateMask(angle.gt(30).And(angle.lt(45)))
end

function get_s1_data(project::String,session,domain::Domain,start_time::String,end_time::String)
    compute_image_url = "https://earthengine.googleapis.com/v1alpha/projects/$project/image:computePixels"
    res = domain.resolution
    xorigin = domain.bbox[1]
    yorigin = domain.bbox[4]
    w,h = domain.shape
    crs = domain.crs
    revisit = 12 # days

    max_retries = 5
    retry_offset = 2 # used to control where to start exponential_backoff
    n_retries = max_retries+retry_offset

    geom = get_ee_region(domain.bbox)

    s1 = (ee.ImageCollection("COPERNICUS/S1_GRD").
        filterDate(start_time,end_time).
        filterBounds(geom).
        filter(ee.Filter.eq("instrumentMode", "IW")).
        filter(ee.Filter.eq("orbitProperties_pass","ASCENDING")).
        filter(ee.Filter.listContains("transmitterReceiverPolarisation", "VV")).
        sort("system:time_start").
        map(s1_qa).
        select("VV")
    )

    dates = get_s1_dates(domain,start_time,end_time)
    n = length(dates)

    jrc = ee.Image("JRC/GSW1_2/GlobalSurfaceWater").select("occurrence").unmask(0).lt(80)

    # n::Int64 = s1.size().getInfo()
    # img_list = s1.toList(n)

    result_arr::Array{Float32} = zeros(Float32,n,h,w,1)
    # result_arr =fill!(Array{Union{Missing, Float64}}(missing,n,h,w,c),-9999)

    @showprogress for i in 1:n
        t1= ee.Date(dates[i])
        t2 = t1.advance(revisit,"day")
        img = ee.Image(s1.filterDate(t1,t2).mean()).
        updateMask(jrc)

        serialized = ee.serializer.encode(img, for_cloud_api=true)

        payload = Dict(
            "expression" => serialized,
            "fileFormat" => "NPY",
            "grid" => Dict(
                "affineTransform" => Dict(
                    "scaleX" => res,
                    "scaleY" => -res,
                    "translateX" => xorigin,
                    "translateY" => yorigin
                ),
                "dimensions" => Dict(
                    "width" => w,
                    "height" => h
                ),
                "crsCode"=> crs
            )
        )

        for j in retry_offset:n_retries
            response = session.post(
                url = compute_image_url,
                data = JSON.json(payload)
            )

            if response.status_code == 200
                arr = response_to_array(response,w,h)
                result_arr[i,:,:,:] = reshape(arr,(h,w,1))
                break
            else
                sleep(exponential_backoff(j))
            end

            if j == n_retries
                println(response.text)
                throw(Base.ProcessFailedException)
                return
            end
        end

    end

    # find no data values and set to missing
    #result_arr[result_arr.==-9999.] .= missing

    return result_arr

end

function get_session(secret_key)
    # import gcloud authentication packages for interacting with the EE REST API
    service_account = pyimport("google.oauth2.service_account")
    gauth = pyimport("google.auth.transport.requests")
    # define the credentials to authorize
    credentials = service_account.Credentials.from_service_account_file(secret_key)
    scoped_credentials = credentials.with_scopes(
        ["https://www.googleapis.com/auth/cloud-platform"])

    # authorize a session to start sending requests
    session = gauth.AuthorizedSession(scoped_credentials)
end

function reof(array::Array{Float32}; variance_threshold::Float64=0.727,n_modes::Int=-1,max_rotations=50)
    spatial_shape = size(array)[2:3]
    flattened = reshape(array,(size(array)[1],prod(spatial_shape)))
    time_shape = size(flattened)[1]

    spatial_mean = mean(flattened,dims=1)
    # time_mean = mean(flattened,dims=2)
    centered = flattened .- spatial_mean

    cov_mat = cov(centered)

    eigen_values,eigen_vectors = eigen(cov_mat)

    # reverse the order of the eigen values and vectors so variance explained is in decreasing order
    eigen_values = eigen_values[end:-1:1]
    eigen_vectors = eigen_vectors[:,end:-1:1]

    if n_modes <= 0
        δλ = eigen_values .* sqrt(2 / time_shape)
        explained_var = cumsum(δλ) / sum(δλ)
        n_modes = sum(explained_var .< variance_threshold)
    end

    eofs = eigen_vectors[:,1:n_modes]

    pcs = Array{Float64}(undef,(time_shape,n_modes))
    for i in 1:n_modes
        pcs[:,i] = centered * eigen_vectors[:,i]
    end

    spatial_modes, temporal_modes = varimax_rotation(eofs,pcs)

    spatial_modes = reshape(spatial_modes,(vcat(spatial_shape...,n_modes)...))

    return spatial_modes, temporal_modes

end

function varimax_rotation(U,A; max_iter=50)
    Z = A * U' # reconstructed data

    #begin varimax rotations
    B = Z * U

    lim = 1
    tot = 1e-6
    cnt = 1

    # need to just initialize V and Ginv
    modes = size(U)[2]
    V = zeros(size(U))
    Ginv = zeros((modes,modes))

    for i in 1:max_iter
    # while lim > tot
        D = diagm(mean(B.^2,dims=1)[1,:])
        C = B.^3
        V = Z' * (C - B * D)
        k,w = eigen((V'*V))
        k = k.^(-1/2)
        for j in 1:size(k)[1]
            if !isfinite(k[j])
                k[j] = 0
            end
        end

        kin = diagm(k)

        Ginv = w * kin * w'

        Uold = U
        U = V * Ginv

        Bold = B
        B = Z*U

        if max((U - Uold)...) < tot
            break
            # return U,B
        end
        cnt += 1
    end

    if cnt >= max_iter
        println("rotation never converged...")
    else
        println("rotation coverged after $cnt iterations")
    end

    U = V * Ginv
    B = Z*U

    return U, B

end

function dualvarimax(Z,Phi)
    U = Phi
    B = Z' * U
    # Unew = zeros(size(Phi))

    lim = 1
    tot = 1e-6
    max_iter = 50
    cnt = 1

    # need to just initialize V and Ginv
    modes = size(U)[2]
    V = zeros(size(U))
    Ginv = zeros((modes,modes))

    for i in 1:max_iter
    # while lim > tot
        D = diagm(mean(B.^2,dims=1)[1,:])
        C = B.^3
        V = Z * (C - B * D)
        k,w = eigen((V'*V))

        kin = diagm(k.^(-1/2))

        Ginv = w * kin * w'

        Uold = U
        U = V * Ginv

        Bold = B
        B = Z'*U

        if max((U - Uold)...) < tot
            break
            # return U,B
        end
        cnt += 1
    end

    U = V* Ginv
    B = Z' * U

    return B, U

end

function get_spt_table(geom=nothing,domain=nothing, use_reach_id=false)
    # define the api service url and endpoint
    spt_api_url = "https://tethys2.byu.edu/localsptapi/api"
    endpoint = "HistoricSimulation"

    if geom === nothing
        geom = get_ee_region(domain.bbox)
    end

    # geom = get_ee_region(domain.bbox)
    reaches = ee.FeatureCollection("users/kelmarkert/public/geoglows/global_drainageline")
    roi_reaches = reaches.filterBounds(geom).sort("to_node",false)
    # need to update filtering process
    # currently gets first reach but needs to be outlet reach
    if use_reach_id
        reach_id = ee.Feature(roi_reaches.first()).get("arcid").getInfo()
        request_url = "$spt_api_url/$endpoint/?reach_id=$reach_id&return_format=json"
    else
        x,y = ee.List(ee.Feature(roi_reaches.first()).geometry().coordinates()).get(0).getInfo()
        request_url = "$spt_api_url/$endpoint/?lon=$x&lat=$y&return_format=json"
    end

    response = HTTP.get(request_url)
    if response.status == 200
        data_dict = JSON.parse(String(response.body))

        df = DataFrame(data_dict["time_series"])
        df[:,"datetime"] = DateTime.(SubString.(df[:,"datetime"],1,19),DateFormat("Y-m-dTH:M:Si"))
        df[:,"flow"] = convert(Array{Float64},df.flow)
    else
        throw(Base.ProcessFailedException)
        return
    end
    return df
end

function match_dates(df::DataFrame,dates::Vector{DateTime})
    # TODO: check to make sure that dates range is within df[:,"datetime"] range
    # produces inconsistent results if not ... not sure why...
    mask = sum(map(x -> x.==df[:,"datetime"],dates),dims=1)[1]
    indices = findall(!iszero,mask)

    df[indices,:]

end

function find_fits(df::DataFrame,sm::Array{Float64},tcp::Array{Float64},stack::Array{Float32})
    # TODO: best fitting???
    # simplisting fitting and will return correlation of tcp to streamflow with rmse of predictions
    x = convert(Array{Float64},df.flow)

    output = Dict{String,Float64}()

    spatial_mean = mean(stack[:,:,:],dims=1)[1,:,:]
    flat_shape = prod(size(spatial_mean))

    errors = Array{Any}(undef,size(df)[1])
    cors = Array{Any}(undef,size(df)[1])

    for i in 1:size(tcp)[2]
        mode_name = "tcp$i"
        y = tcp[:,i]

        for j in 1:3
            model = fit(x,y,j)

            test_preds = synthesize(model,df.flow,sm[:,:,i],spatial_mean)

            for k in 1:size(test_preds)[1]
                diff = test_preds[k,:,:] - stack[k,:,:,1]
                errors = mean(abs.(diff))
                cors = cor(reshape(test_preds[k,:,:],flat_shape),reshape(stack[k,:,:,1],flat_shape))
            end

            push!(output, "$(mode_name)_fit_correlation_order$j" => cor(y,model.(x)))
            # push!(output, "$(mode_name)_pred_correlation_order$j" => mean(cors))
            push!(output, "$(mode_name)_pred_rmse_order$j"=>mean(errors))
        end

    end

    return output

end

function synthesize(model,streamflow,spatial_mode,spatial_mean)
    z = model.(streamflow)

    preds = Array{Float64}(undef,(vcat(length(z),size(spatial_mode)...)...))

    for i in 1:length(z)
        preds[i,:,:] = z[i]*spatial_mode+spatial_mean
    end

    return preds
end

function persist_data(fname::String, dataset, domain, dates)
    # TODO: write data persistance function to save either sar or synthesized data

end

# end module
end


text:
"get mole fraction (X) from mass fraction (Y)"
function Y2X(gas, Y)
    return Y ./ gas.MW / sum(Y./gas.MW)
end
Y2X(gas, Y, mean_MW) = Y2X(gas, Y)
export Y2X

"get concentration (C) from mass fraction (Y)"
function Y2C(gas, Y, ρ_mass)
    return Y * ρ_mass ./ gas.MW
end
export Y2C

"get mole fraction (X) from concentration (C)"
function C2X(gas, C)
    return C ./ sum(C)
end
export C2X

"get concentration (C) from mole fraction (X)"
function X2C(gas, X, ρ_mass)
    return X * ρ_mass / sum(X.*gas.MW)
end
export X2C

"get mass fraction (Y) from mole fraction (X)"
function X2Y(gas, X)
    return X .* gas.MW / sum(X.*gas.MW)
end
X2Y(gas, X, mean_MW) = X2Y(gas, X)
export X2Y

"""
Get species index of a species
# Example
```
    species_index(gas, "O2")
```
"""
function species_index(gas, species)
    return findfirst(gas.species_names .== species)
end
export species_index

"""
"mainly for testing code, will be removed in the future"
"please customize such functions following this example"
"""
function set_states(gas::Solution, T, P, Y)
    mean_MW = 1.0 / dot(Y, 1 ./ gas.MW)
    ρ_mass = P / R / T * mean_MW
    X = Y2X(gas, Y, mean_MW)
    C = Y2C(gas, Y, ρ_mass)
    cp_mole, cp_mass = get_cp(gas, T, X, mean_MW)
    h_mole = get_H(gas, T, Y, X)
    S0 = get_S(gas, T, P, X)
    wdot = wdot_func(gas.reaction, T, C, S0, h_mole)
    return wdot
end
export set_states


text:
# 2D FWI gradient test with 4 sources
# The receiver positions and the source wavelets are the same for each of the four experiments.
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#

using PyCall, PyPlot, opesciSLIM.TimeModeling

## Set up model structure
n = (80,80)	# (x,y,z) or (x,z)
d = (10.,10.)
o = (0.,0.)

# Velocity [km/s]
v = ones(Float32,n) + 0.4f0
v[:,Int(round(end/2)):end] = 2.5f0
v0 = smooth10(v,n)

# Slowness squared [s^2/km^2]
m = (1f0./v).^2
m0 = (1f0./v0).^2
dm = m0 - m

# Setup info and model structure
nsrc = 1	# number of sources
ntComp = 250
info = Info(prod(n), nsrc, ntComp)	# number of gridpoints, number of experiments, number of computational time steps
model = Model(n,d,o,m)
model0 = Model(n,d,o,m0)

## Set up receiver geometry
nxrec = 101
xrec = linspace(50f0,750f0,nxrec)
yrec = 0.
zrec = linspace(40f0,40f0,nxrec)

# receiver sampling and recording time
timeR = 1200f0	# receiver recording time [ms]
dtR = calculate_dt(n,d,o,v)	# receiver sampling interval

# Set up receiver structure
recGeometry = Geometry(xrec,yrec,zrec;dt=dtR,t=timeR,nsrc=nsrc)

## Set up source geometry (cell array with source locations for each shot)
xsrc = convertToCell([400f0])
ysrc = convertToCell([0f0])
zsrc = convertToCell([20f0])

# source sampling and number of time steps
timeS = 1200f0
dtS = calculate_dt(n,d,o,v)	# receiver sampling interval

# Set up source structure
srcGeometry = Geometry(xsrc,ysrc,zsrc;dt=dtS,t=timeS)

# setup wavelet
f0 = 0.01f0
wavelet = source(timeS,dtS,f0)

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

# Gradient test
h = 1.f-1
iter = 10
error1 = zeros(iter)
error2 = zeros(iter)
h_all = zeros(iter)
srcnum = 1:nsrc
modelH = deepcopy(model0)

# Observed data
F = joModeling(info,model,srcGeometry,recGeometry)
q = joData(srcGeometry,wavelet)
d = F*q

# FWI gradient and function value for m0
Jm0, grad = fwi_objective(model0,q,d)

for j=1:iter
	# FWI gradient and function falue for m0 + h*dm
	modelH.m = model0.m + h*dm
	Jm, gradm = fwi_objective(modelH,q,d)

	dJ = dot(grad,vec(dm))

	# Check convergence
	error1[j] = abs(Jm - Jm0)
	error2[j] = abs(Jm - (Jm0 + h*dJ))

	println(h, " ", error1[j], " ", error2[j])
	h_all[j] = h
	h = h/2f0
end

# Plot errors
loglog(h_all, error1); loglog(h_all, 1e2*h_all)
loglog(h_all, error2); loglog(h_all, 1e2*h_all.^2)
legend([L"$\Phi(m) - \Phi(m0)$", "1st order", L"$\Phi(m) - \Phi(m0) - \nabla \Phi \delta m$", "2nd order"], loc="lower right")
xlabel("h")
ylabel(L"Error $||\cdot||^\infty$")
title("FWI gradient test")
#axis((h_all[end], h_all[1], 1.0e-8,500))

max_stars_repo_name: SINBADconsortium/opesciSLIM

max_issues_repo_name: SINBADconsortium/opesciSLIM

max_forks_repo_name: SINBADconsortium/opesciSLIM


text:
@cache mutable struct ISSEMCache{uType,rateType,N,noiseRateType,randType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  gtmp::noiseRateType
  gtmp2::rateType
  nlsolver::N
  dW_cache::randType
end

function alg_cache(alg::ISSEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,
                   ::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  γ, c = alg.theta,zero(t)
  nlsolver = OrdinaryDiffEq.build_nlsolver(alg,u,uprev,p,t,dt,f,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,γ,c,Val(true))
  fsalfirst = zero(rate_prototype)
  gtmp = zero(noise_rate_prototype)
  if is_diagonal_noise(prob)
    gtmp2 = gtmp
    dW_cache = nothing
  else
    gtmp2 = zero(rate_prototype)
    dW_cache = zero(ΔW)
  end

  ISSEMCache(u,uprev,fsalfirst,gtmp,gtmp2,nlsolver,dW_cache)
end

mutable struct ISSEMConstantCache{N} <: StochasticDiffEqConstantCache
  nlsolver::N
end

function alg_cache(alg::ISSEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,
                   ::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  γ, c = alg.theta,zero(t)
  nlsolver = OrdinaryDiffEq.build_nlsolver(alg,u,uprev,p,t,dt,f,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,γ,c,Val(false))
  ISSEMConstantCache(nlsolver)
end

@cache mutable struct ISSEulerHeunCache{uType,rateType,N,noiseRateType,randType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  gtmp::noiseRateType
  gtmp2::rateType
  gtmp3::noiseRateType
  nlsolver::N
  dW_cache::randType
end

function alg_cache(alg::ISSEulerHeun,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,
                   ::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  γ, c = alg.theta,zero(t)
  nlsolver = OrdinaryDiffEq.build_nlsolver(alg,u,uprev,p,t,dt,f,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,γ,c,Val(true))
  fsalfirst = zero(rate_prototype)

  gtmp = zero(noise_rate_prototype)
  gtmp2 = zero(rate_prototype)

  if is_diagonal_noise(prob)
      gtmp3 = gtmp2
      dW_cache = nothing
  else
      gtmp3 = zero(noise_rate_prototype)
      dW_cache = zero(ΔW)
  end

  ISSEulerHeunCache(u,uprev,fsalfirst,gtmp,gtmp2,gtmp3,nlsolver,dW_cache)
end

mutable struct ISSEulerHeunConstantCache{N} <: StochasticDiffEqConstantCache
  nlsolver::N
end

function alg_cache(alg::ISSEulerHeun,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,
                   ::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  γ, c = alg.theta,zero(t)
  nlsolver = OrdinaryDiffEq.build_nlsolver(alg,u,uprev,p,t,dt,f,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,γ,c,Val(false))
  ISSEulerHeunConstantCache(nlsolver)
end

max_stars_repo_name: JuliaDiffEq/StochasticDiffEq


text:
module HypergraphsEpidemics

using CSV
using DataFrames
using Query
using Tables
using Dates
using DataStructures
using Distributions
using SimpleHypergraphs
using Statistics
using Combinatorics
using Pipe: @pipe
using Random
using JSON
using JSON3
using Serialization

export Person
export Household
export Company
export Contact_simulation_options

export generate_dataset
export start_simulation

include("common/Constants.jl")
include("common/EventWriter.jl")
include("common/utils.jl")
include("common/structs.jl")
include("RoutineAutomaton.jl")
include("DatasetGeneration.jl")
include("Simulation.jl")

include("../deps/HGEpidemics/src/HGEpidemics.jl")

using .HGEpidemics

end

max_stars_repo_name: Luigi-Crisci/Hypergraphs-Epidemics

max_issues_repo_name: Luigi-Crisci/Hypergraphs-Epidemics

max_forks_repo_name: Luigi-Crisci/Hypergraphs-Epidemics


text:
#!/usr/bin/env julia

# script with a simple command line interface

# the first line allows the script to be run like an executable from the
# terminal

# The "ArgParse" framework provies a simpl framework for parsing and type
# checking arguments passed from the terminal
using ArgParse

function main(args::Dict{String,Any})
    println("string: $(args["string"])")
    println("number: $(args["number"])")
    println("flag:   $(args["flag"])")
end

# this isinteractive() check makes sure that main is not run when loading
# this script into the REPL
if !isinteractive()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--string", "-s"
            help = "some string data"
            required = true
        "--number", "-n"
            help = "any integer will do"
            arg_type = Int
            default = 7
        "--flag", "-f"
            help = "a boolean flag"
            action = :store_true
    end

    # parses julia's ARGS list based on the arg table data
    main(parse_args(s))
end

# test with terminal commands like
#
#   ./cli-argparse.jl
#   ./cli-argparse.jl -h
#   ./cli-argparse.jl -s foo
#

# test with REPL commands like
#
#   include("cli-argparse.jl")
#   main(Dict("string" => "foo", "number" => 4, "flag" => false))
#

max_stars_repo_name: ccoffrin/julia-examples

max_issues_repo_name: ccoffrin/julia-examples

max_forks_repo_name: ccoffrin/julia-examples


text:
# Model
include_model("walker")

# Visualize
# - Pkg.add any external deps from visualize.jl
include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)

"""
Useful methods

joint-space inertia matrix:
	M_func(model, q)

position (px, pz) of foot 1 heel:
	kinematics_2(model,
		get_q⁺(x), body = :calf_1, mode = :ee)

position (px, pz) of foot 2 heel:
	kinematics_2(model,
		get_q⁺(x), body = :calf_2, mode = :ee)

position (px, pz) of foot 1 toe:
		kinematics_3(model,
	        get_q⁺(x), body = :foot_1, mode = :ee)

position (px, pz) of foot 2 toe:
		kinematics_3(model,
	        get_q⁺(x), body = :foot_2, mode = :ee)

jacobian to foot 1 heel
	jacobian_2(model, get_q⁺(x), body = :calf_1, mode = :ee)

jacobian to foot 2 heel
	jacobian_2(model, get_q⁺(x), body = :calf_2, mode = :ee)

jacobian to foot 1 toe
	jacobian_3(model, get_q⁺(x), body = :foot_1, mode = :ee)

jacobian to foot 2 toe
	jacobian_3(model, get_q⁺(x), body = :foot_2, mode = :ee)

gravity compensating torques
	initial_torque(model, q, h)[model.idx_u] # returns torques for system
	initial_torque(model, q, h) # returns torques for system and contact forces

	-NOTE: this method is iterative and is not differentiable (yet)
		-TODO: implicit function theorem to compute derivatives

contact models
	- maximum-dissipation principle:
		contact_constraints(model, T)
	- no slip:
		contact_no_slip_constraints(T)

loop constraint
	- constrain variables at two indices to be the same:
		loop_constraints(model, x_idx, t1_idx, t2_idx)
"""

function get_q⁺(x)
	view(x, model.nq .+ (1:model.nq))
end

# Horizon
T = 26
Tm = convert(Int, floor(0.5 * T))

# Time step
tf = 2.5
h = tf / (T - 1)

# Configurations
# 1: x pos
# 2: z pos
# 3: torso angle (rel. to downward vertical)
# 4: thigh 1 angle (rel. to downward vertical)
# 5: calf 1 (rel. to downward vertical)
# 6: thigh 2 (rel. to downward vertical)
# 7: calf 2 (rel. to downward vertical)
# 8: foot 1 (rel. to downward vertical)
# 9: foot 2 (rel. to downward vertical)
function initial_configuration(model, θ_torso, θ_thigh_1, θ_leg_1, θ_thigh_2)
    q1 = zeros(model.nq)
    q1[3] = θ_torso
    q1[4] = θ_thigh_1
    q1[5] = θ_leg_1
    z1 = model.l_thigh1 * cos(q1[4]) + model.l_calf1 * cos(q1[5])
    # q1[6] = θ_thigh_2
    # q1[7] = -1.0 * acos((z1 - model.l_thigh2 * cos(q1[6])) / model.l_calf2)
	q1[6] = θ_thigh_1
	q1[7] = θ_leg_1
	q1[2] = z1
	q1[8] = pi / 2.0
	q1[9] = pi / 2.0

    p1 = kinematics_2(model, q1, body = :calf_1, mode = :ee)
    p2 = kinematics_2(model, q1, body = :calf_2, mode = :ee)
    @show stride = abs(p1[1] - p2[1])

    q1[1] = -0.75#-1.0 * p1[1]

    qM = copy(q1)
    qM[4] = q1[6]
    qM[5] = q1[7]
    qM[6] = q1[4]
    qM[7] = q1[5]
    qM[1] = abs(p2[1])

    pM_1 = kinematics_2(model, qM, body = :calf_1, mode = :ee)
    pM_2 = kinematics_2(model, qM, body = :calf_2, mode = :ee)

    qT = copy(q1)
    qT[1] = 0.75#q1[1] + 1.0#2 * stride

    pT_1 = kinematics_2(model, qT, body = :calf_1, mode = :ee)
    pT_2 = kinematics_2(model, qT, body = :calf_2, mode = :ee)

    return q1, qM, qT
end

q1, qM, qT = initial_configuration(model, -pi / 100.0, pi / 20.0, -pi / 10.0, -pi / 25.0)
q_ref = linear_interpolation(q1, qT, T)
x_ref = configuration_to_state(q_ref)
visualize!(vis, model, q_ref, Δt = h)

# Control
# u = (τ1..7, λ1..4, β1..8, ψ1..4, η1...8, s1)
# τ1: torso angle
# τ2: thigh 1 angle
# τ3: calf 1
# τ4: thigh 2
# τ5: calf 2
# τ6: foot 1
# τ7: foot 2

# ul <= u <= uu
u1 = initial_torque(model, q1, h)[model.idx_u] # gravity compensation for current q
_uu = Inf * ones(model.m)
_uu[model.idx_u] .= 10.0
_ul = zeros(model.m)
_ul[model.idx_u] .= -10.0
ul, uu = control_bounds(model, T, _ul, _uu)

qL = [-Inf; -Inf; q1[3:end] .- pi / 4.0; -Inf; -Inf; q1[3:end] .- pi / 4.0]
qU = [Inf; q1[2] + 0.001; q1[3:end] .+ pi / 4.0; Inf; Inf; q1[3:end] .+ pi / 4.0]

q2 = copy(q1)
q2[1] = q1[1] + 0.25 * h
xl, xu = state_bounds(model, T,
    qL, qU,
    x1 = [q1; q2], # initial state
    xT = [Inf * ones(model.nq); qT]) # goal state

# Objective
include_objective(["velocity", "nonlinear_stage", "control_velocity"])

x0 = configuration_to_state(q_ref)

# penalty on slack variable
obj_penalty = PenaltyObjective(1.0e4, model.m)

# quadratic tracking objective
# Σ (x - xref)' Q (x - x_ref) + (u - u_ref)' R (u - u_ref)
q_penalty = 0.1 * ones(model.nq)
q_penalty[1] = 0.0
q_penalty[2] = 10.0
q_penalty[3] = 1.0
x_penalty = [q_penalty; q_penalty]
obj_control = quadratic_tracking_objective(
    [h * Diagonal(x_penalty) for t = 1:T],
    [h * Diagonal([1.0 * ones(model.nu)..., 0.0 * ones(model.m - model.nu)...]) for t = 1:T-1],
    [x_ref[T] for t = 1:T],
    [[zeros(model.nu); zeros(model.m - model.nu)] for t = 1:T-1])

# quadratic velocity penalty
# Σ v' Q v
q_v = 10.0 * ones(model.nq)
# q_v[1] = 0.0
# q_v[2] = 0.0
obj_velocity = velocity_objective(
    [h * Diagonal(q_v) for t = 1:T-1],
    model.nq,
    h = h)#,
    # idx_angle = collect([3, 4, 5, 6, 7, 8 ,9]))

function l_foot_height(x, u, t)
	ϕ_des = [0.1; 0.1; 0.1; 0.1]
	ϕ = ϕ_func(model, get_q⁺(x))

	idx = (ϕ .< ϕ_des)

	return 1000.0 * (ϕ[idx] - ϕ_des[idx])' * (ϕ[idx] - ϕ_des[idx])
end

l_foot_height(x) = l_foot_height(x, nothing, T)

obj_foot_height = nonlinear_stage_objective(l_foot_height, l_foot_height)

obj = MultiObjective([obj_penalty,
                      obj_control,
                      obj_velocity,
					  obj_foot_height])
					  # obj_control_velocity])
                      # obj_th])#,
                      # obj_fh1,
                      # obj_fh2,
					  # obj_tf])
					  # obj_forward])

# Constraints
include_constraints(["contact"])#, "loop", "free_time"])
# con_loop = loop_constraints(model, collect([(2:7)...,(9:14)...]), 1, T)
con_contact = contact_constraints(model, T)
# con_free_time = free_time_constraints(T)
con = multiple_constraints([con_contact])#, con_free_time, con_loop])

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# trajectory initialization
u0 = [0.01 * randn(model.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

# Solve
# NOTE: run multiple times to get good trajectory
include_snopt()
@time z̄, info = solve(prob, copy(z0),
    nlp = :SNOPT7,
    tol = 1.0e-4, c_tol = 1.0e-3, mapl = 5,
    time_limit = 60 * 1)

@show check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)
vis = Visualizer()
render(vis)
visualize!(vis, model,
	[[x̄[1][1:model.nq] for i = 1:10]...,
	 state_to_configuration(x̄)...,
	 [x̄[T][model.nq .+ (1:model.nq)] for i = 1:10]...], Δt = h)
open(vis)
# @save joinpath(pwd(), "examples/trajectories/walker_steps.jld2") z̄

# using Plots
# fh1 = [kinematics_3(model,
#     state_to_configuration(x̄)[t], body = :foot_1, mode = :com)[2] for t = 1:T]
# fh2 = [kinematics_3(model,
#     state_to_configuration(x̄)[t], body = :foot_2, mode = :com)[2] for t = 1:T]
# plot(fh1, linetype = :steppost, label = "foot 1")
# plot!(fh2, linetype = :steppost, label = "foot 2")

# plot(hcat(ū...)[1:7, :]',
#     linetype = :steppost,
#     label = "",
#     color = :red,
#     width = 2.0)
#
# plot!(hcat(u_proj...)[1:4, :]',
#     linetype = :steppost,
#     label = "",
#     color = :black)
#
# plot(hcat(u_proj...)[5:6, :]',
#     linetype = :steppost,
#     label = "",
#     width = 2.0)

# plot(hcat(state_to_configuration(x̄)...)'[1:3],
#     color = :red,
#     width = 2.0,
#     label = "")
#
# plot!(hcat(state_to_configuration(x_proj)...)',
#     color = :black,
#     label = "")


text:
@testset "kaggle dataset" begin
    dataset = "ningjingyu/fetchtest"

    f = kdownload(dataset, pwd())
    DataDeps.unpack(f)

    open(joinpath(pwd(), "FetchTest", "FetchTest.txt"), "r") do file
        @test readline(file) == "Test"
    end

    rm(joinpath(pwd(), "FetchTest"), recursive=true, force=true)
end

@testset "kaggle url" begin
    urls = [
        "https://www.kaggle.com/ningjingyu/fetchtest",
        "https://www.kaggle.com/ningjingyu/fetchtest/tasks",
        "https://www.kaggle.com/ningjingyu/fetchtest/code",
        "https://www.kaggle.com/ningjingyu/fetchtest/discussion",
        "https://www.kaggle.com/ningjingyu/fetchtest/activity",
        "https://www.kaggle.com/ningjingyu/fetchtest/metadata",
        "https://www.kaggle.com/ningjingyu/fetchtest/settings",
    ]

    for url in urls
        f = kdownload(url, pwd())
        DataDeps.unpack(f)

        open(joinpath(pwd(), "FetchTest", "FetchTest.txt"), "r") do file
            @test readline(file) == "Test"
        end

        rm(joinpath(pwd(), "FetchTest"), recursive=true, force=true)
    end
end


text:
# Dataset types and utility functions to construct datasets

searchdir(path, key) = filter(x -> contains(x, key), readdir(path))

# TODO(ryuichi): should be more generic
function push_delta(src)
    D, T = size(src)
    src = repmat(src, 2)
    for t=2:T-1
        src[D+1:end,t] = -0.5*src[1:D,t-1] + 0.5*src[1:D,t+1]
    end
    src
end

abstract Dataset

immutable ParallelDataset <: Dataset
    X::Matrix{Float64}
    Y::Matrix{Float64}

    # keep statistics in case of recovering original data from
    # standarized data
    Xmean::Matrix{Float64}
    Xstd::Matrix{Float64}
    Ymean::Matrix{Float64}
    Ystd::Matrix{Float64}

    diff::Bool

    function ParallelDataset(path;
                             diff::Bool=false,
                             joint::Bool=true,
                             standarize::Bool=false,
                             ignore0th::Bool=true,
                             add_delta::Bool=false,
                             suffix::AbstractString="_parallel.jld",
                             keepstat::Bool=false,
                             nmax::Int=100)
        files = searchdir(path, suffix)
        sort!(files)

        println("$(length(files)) training data found")

        XY = ones(1,1)
        totalframes::Int = 0
        totalphrases::Int = 0

        count = 0
        for filename in files
            # TODO(ryuichi) allow costom file format?
            f = load(joinpath(path, filename))
            src, tgt = f["src"], f["tgt"]
            src_x, tgt_x = src["feature_matrix"], tgt["feature_matrix"]

            if ignore0th
                src_x, tgt_x = src_x[2:end,:], tgt_x[2:end,:]
            end

            if add_delta
                src_x = push_delta(src_x)
                tgt_x = push_delta(tgt_x)
            end

            # use differencial
            if diff
                tgt_x = tgt_x - src_x
            end

            # Create joint features matries of source and target speaker
            combined = vcat(src_x, tgt_x)
            @assert size(combined) == (size(src_x, 1)*2, size(src_x, 2))

            if totalframes == 0
                XY = combined
            else
                XY = hcat(XY, combined)
            end

            totalframes += size(combined, 2)
            totalphrases += 1

            count += 1
            count >= nmax && break
        end

        println("total number of frames: $(totalframes)")
        println("total number of phrases: $(totalphrases)")

        X = ones(1,1)
        Y = ones(1,1)
        if joint
            X = XY
            Y = ones(1,1) # not used
        else
            D = div(size(XY, 1), 2)
            X = XY[1:D,:]
            Y = XY[D+1:end,:]
        end

        Xmean = mean(X,2)
        Xstd = std(X,2)
        Ymean = mean(Y,2)
        Ystd = std(Y,2)
        if standarize
            X = (X - Xmean) / Xstd
            Y = (Y - Ymean) / Ystd
        end

        # avoid nan
        if joint
            Ymean = ones(1,1)
            Ystd = ones(1,1)
        end

        @assert all(isfinite.(X))
        @assert all(isfinite.(Y))
        @assert all(isfinite.(Xmean))
        @assert all(isfinite.(Xstd))
        @assert all(isfinite.(Ymean))
        @assert all(isfinite.(Ystd))

        if !keepstat
            Xmean = ones(1,1)
            Xstd = ones(1,1)
            Ymean = ones(1,1)
            Ystd = ones(1,1)
        end

        new(X, Y, Xmean, Xstd, Ymean, Ystd, diff)
    end
end

# GVDataset represents a global variance dataset.
immutable GVDataset <: Dataset
    X::Matrix{Float64}

    function GVDataset(path;
                       ignore0th::Bool=true,
                       add_delta::Bool=false,
                       suffix::AbstractString=".jld",
                       nmax::Int=100)
        files = searchdir(path, suffix)
        sort!(files)

        println("$(length(files)) training data found")

        X = zeros(0, 0)
        totalphrases::Int = 0

        for filename in files
            # TODO(ryuichi) allow costom file format?
            f = load(joinpath(path, filename))
            tgt = f["feature_matrix"]

            if ignore0th
                tgt = tgt[2:end,:]
            end

            if add_delta
                tgt = push_delta(tgt)
            end

            gv = var(tgt, 2)
            if sum(isnan.(gv)) == 0
                if totalphrases == 0
                    X = gv
                else
                    X = hcat(X, gv)
                end
            end

            totalphrases += 1

            totalphrases >= nmax && break
        end

        println("total number of phrases: $(totalphrases)")

        @assert all(isfinite.(X))

        new(X)
    end
end


text:
module RiskyCommunity

    # __precompile__(true)
    using Parameters
    using JuMP
    using AxisArrays

    include("base_structs.jl")
    include("community_agents.jl")
    include("DER_assets.jl")
    include("grid_agents.jl")
    include("opt_models.jl")

    export build_model, build_test_model, build_VAR_model, add_objective, solve_model, fix_model_variable!, read_results
    export Period, Scenario, Result, Demand, Price
    export DER_Agent, Battery
    export Battery, PV
    export Conv_Generator, Res_Generator, Load_Shed, Central_Opt, Aggregator

end # module


text:
#=##############################################################################
# DESCRIPTION
    Functions for exporting geometries into PLOT3D xyz format
# AUTHORSHIP
  * Author    : Eduardo J Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Jul 2019
  * License   : MIT License
=###############################################################################

"""
Generate a structured PLOT3D mesh
"""
function generatePLOT3D(grid::GridExtentions, name::String; path="",
                        num=nothing, ext="xyz", trickwopwop=false)

    fname = joinpath(path, name*(num!=nothing ? ".$num" : "")*"."*ext)
    f = open(fname, "w")

    ndivs = get_ndivsnodes(grid)
    dims = grid.dims

    # imax jmax kmax
    if trickwopwop
        print(f, ndivs[1], " ", dims==3  ? ndivs[2]*ndivs[3] :  ndivs[2], " 1")
    else
        for (i,ndiv) in enumerate(ndivs)
            print(f, i!=1 ? " " : "", ndiv)
        end
        for i in 1:3-dims
            print(f, i!=1 ? " " : "", 1)
        end
    end
    print(f, "\n")

    # imax × jmax ×kmax floating point x coordinates
    for n in 1:grid.nnodes
        s = @sprintf "%8.6f" get_node(grid, n)[1]
        print(f, n!=1 ? (n%10==0 ? "\n" : " ") : "", s)
    end
    print(f, "\n")
    # imax × jmax ×kmax floating point y coordinates
    for n in 1:grid.nnodes
        s = @sprintf "%8.6f" get_node(grid, n)[2]
        print(f, n!=1 ? (n%10==0 ? "\n" : " ") : "", s)
    end
    print(f, "\n")
    # imax × jmax ×kmax floating point z coordinates
    for n in 1:grid.nnodes
        s = @sprintf "%8.6f" get_node(grid, n)[3]
        print(f, n!=1 ? (n%10==0 ? "\n" : " ") : "", s)
    end
    close(f)
end


text:
## -------------------------------------------------------------------
# model vs exp comparison
let
    METHODS = [
        :FBA_Z_FIX_MAX_VG_MIN_COST,
        :FBA_Z_FIX_MAX_VG_MAX_COST,
        :FBA_Z_FIX_MIN_VG_MIN_COST,
        # :FBA_Z_FIX_MIN_VG_MAX_COST,
        # :FBA_Z_FIX_ADJ_MAX_VG_MIN_COST,
        # :FBA_Z_FIX_ADJ_MAX_VG_MAX_COST,
        # :FBA_Z_FIX_ADJ_MIN_VG_MIN_COST,
        # :FBA_Z_FIX_ADJ_MIN_VG_MAX_COST,
        # :FBA_Z_FIX_MAX_ATP_MIN_COST,
        :FBA_Z_FIX_MAX_ATP_MAX_COST,
        :ME_MAX_POL,
        # :ME_MAX_POL_B0
    ]

    colors = Dict(
        :exp => :black,
        :FBA_Z_FIX_MAX_VG_MIN_COST => :blue,
        :FBA_Z_FIX_MAX_VG_MAX_COST => :red,
        :FBA_Z_FIX_MIN_VG_MIN_COST => :black,
        :FBA_Z_FIX_MIN_VG_MAX_COST => :red,
        :FBA_Z_FIX_MAX_ATP_MIN_COST => :blue,
        :FBA_Z_FIX_MAX_ATP_MAX_COST => :black,
        :ME_MAX_POL => :yellow,
        :ME_MAX_POL_B0 => :purple,
    )

    ms = 8

    # [:none, :auto, :circle, :rect, :star5, :diamond, :hexagon, :cross, :xcross, :utriangle,
    # :dtriangle, :rtriangle, :ltriangle, :pentagon, :heptagon, :octagon, :star4, :star6,
    # :star7, :star8, :vline, :hline, :+, :x].
    shapes = Dict(
        :exp => :star,
        :FBA_Z_FIX_MAX_VG_MIN_COST => :utriangle,
        :FBA_Z_FIX_MAX_VG_MAX_COST => :utriangle,
        :FBA_Z_FIX_MIN_VG_MIN_COST => :square,
        :FBA_Z_FIX_MIN_VG_MAX_COST => :square,
        :FBA_Z_FIX_MAX_ATP_MIN_COST => :hexagon,
        :FBA_Z_FIX_MAX_ATP_MAX_COST => :hexagon,
        :ME_MAX_POL => :circle,
        :ME_MAX_POL_B0 => :square,
    )

    iJR = ChN.iJR904
    Data = ChN.NanchenData

    Nd_rxns_map = iJR.load_rxns_map2() # inner reacts
    iJR_ider_subs = iJR.load_inner_rxns_subs()

    sort_iders(iders) = sort!(deepcopy(iders),
        by = (ider) -> ider == "GLC" ? "AAA" : ider
    )

    # Collect
    src = nameof(Data)

    DAT = ChE.load_DAT(src)

    EXCH_FLX_IDERS = sort_iders(DAT[:EXCH_FLX_IDERS])
    # EXPS = DAT[:EXPS]

    # experiments with approx D
    Dintervals = Dict(
        0.05 => 1:3,
        0.1 => 4:9,
        0.15 => 4:11,
        0.2 => 10:11,
        0.25 => 10:15,
        0.3 => 12:15,
        0.35 => 12:17,
        0.4 => 16:17,
        :tot => 1:17
    )

    @info("Doing", src)

    for (avD, exps) in Dintervals
        dat_pool = Dict()
        # Experiments
        edat = get!(dat_pool, :exp, Dict())
        for exp in exps
            !haskey(DAT, :exp, :flx, "GLC", exp) && continue
            glc_exp_val = DAT[:exp, :flx, "GLC", exp]

            for (sub, iJR_iders) in iJR_ider_subs
                for iJR_ider in iJR_iders
                    Nd_ider = Nd_rxns_map[iJR_ider]
                    exp_val = abs(DAT[:exp, :flx, Nd_ider, exp] / glc_exp_val)

                    # @info(:exp, exp, iJR_ider, exp_val); println()

                    eidat = get!(edat, iJR_ider, [])
                    push!(eidat, exp_val)
                end
            end
        end

        # methods
        for method in METHODS
            mdat = get!(dat_pool, method, Dict())

            for exp in exps
                !haskey(DAT, :exp, :flx, "GLC", exp) && continue
                glc_exp_val = DAT[:exp, :flx, "GLC", exp]

                @info(method, exp); println()

                for (sub, iJR_iders) in iJR_ider_subs
                    (sub in ["others"]) && continue
                    for iJR_ider in iJR_iders
                        !haskey(DAT, method, :flx, iJR_ider, exp) && continue
                        model_val = abs(DAT[method, :flx, iJR_ider, exp] / glc_exp_val)

                        @info(method, exp, iJR_ider, model_val); println()

                        midat = get!(mdat, iJR_ider, [])
                        push!(midat, model_val)
                    end
                end
            end
        end

        # Plots
        xticks_pos = []
        xticks_iders = []
        fontsize = 13
        p = plot(;xlabel = "Reactions",
            ylabel = "flux / glc exch.", xgrid = false,
            # dpi = 1000,
            thickness_scaling = 1.3,
            xguidefontsize = fontsize, yguidefontsize = fontsize
        )
        sep1 = 5
        sep2 = 3 * sep1

        # backgound
        gc = 1
        last_gc = gc
        ci = 1
        for sub in NiJR.subs()
            (sub in ["others"]) && continue
            iJR_iders = iJR_ider_subs[sub]

            for iJR_ider in iJR_iders
                vline!(p, [gc - sep1];
                    label = "", color = :darkgray, alpha = 0.8, ls = :dash
                )
                gc += sep1
                for method in METHODS
                    gc += sep1
                end
                gc += sep2
            end
            fillalpha = isodd(ci) ? 0.6 : 0.3
            vspan!(p, [last_gc - sep1, gc - sep1]; label = "",
                linecolor = :darkgray, fillcolor = :darkgray,
                fillalpha, linealpha = 0.6
            )
            last_gc = gc
            ci += 1
        end

        # fluxes
        gc = 1
        # f(v) = iszero(v) ? -6.0 : log(v)
        f(v) = v
        for sub in NiJR.subs()
            (sub in ["others"]) && continue

            iJR_iders = iJR_ider_subs[sub]
            for iJR_ider in iJR_iders

                # Experimental
                edat = dat_pool[:exp]
                !haskey(edat, iJR_ider) && continue
                m = mean(abs.(edat[iJR_ider]))
                e = std(abs.(edat[iJR_ider]))

                @info(iJR_ider, m, e); println()

                scatter!(p, [gc + sep1], [f(m)];
                    yerr = [e],
                    color = colors[:exp], label = "",
                    m = (ms, shapes[:exp]),
                    alpha = 0.8
                )
                push!(xticks_pos, gc + 2 * sep1)
                push!(xticks_iders, iJR_ider)
                gc += sep1

                # modeled
                for method in METHODS
                    mdat = dat_pool[method]
                    color = colors[method]
                    μ = mean(abs.(mdat[iJR_ider]))
                    e = std(abs.(mdat[iJR_ider]))
                    m = (ms, get(shapes, method, :square))
                    scatter!(p, [gc + sep1], [f(μ)];
                        yerr = [e],
                        color, label = "", m, alpha = 0.8
                    )
                    gc += sep1
                end
                gc += sep2
            end
        end

        xticks = (xticks_pos, xticks_iders)
        xrotation = 35
        size = (1100, 500)
        thickness_scaling = 1.3
        plot!(p; xticks, xrotation, size, thickness_scaling)
        mysavefig(p, "inner_flx_comparizon"; avD)
        PS[:inner_flx_comparizon, :all_subs, avD] = deepcopy(p)
    end
end


text:
"""
    IntSet(n)

Subset of `{1,...,n}` with `O(1)` `in()` and `push!()`.
"""
struct IntSet
    set::Vector{Int}
    iselement::Vector{Bool}
end

IntSet(n) = IntSet(sizehint!(Int[],n), fill(false,n))

limit(s::IntSet) = length(s.iselement)
Base.length(s::IntSet) = length(s.set)
Base.iterate(s::IntSet, args...) = iterate(s.set, args...)
Base.@propagate_inbounds Base.in(v::Int, s::IntSet) = s.iselement[v]

function Base.push!(s::IntSet, v::Int)
    @boundscheck @assert 1 <= v <= limit(s)
    @inbounds begin
        if !s.iselement[v]
            push!(s.set,v)
            s.iselement[v] = true
        end
        return s
    end
end

Base.@propagate_inbounds function Base.union!(s::IntSet, itr)
    for v in itr
        push!(s,v)
    end
    return s
end

function Base.empty!(s::IntSet)
    @inbounds @. s.iselement[s.set] = false
    empty!(s.set)
    return s
end

Base.sort!(s::IntSet) = sort!(s.set)

max_stars_repo_name: ettersi/SparseFactorizations

max_issues_repo_name: ettersi/SparseFactorizations

max_forks_repo_name: ettersi/SparseFactorizations


text:
using Test, VAN

@testset "test product model" begin
    include("psa.jl")
end

@testset "test autoregressive model" begin
    include("arm.jl")
end


text:
import HTTP, CSV, JSON
using Dashboards, PlotlyJS, DataFrames

external_stylesheets = ["https://codepen.io/chriddyp/pen/bWLwgP.css"]
println("downloading data for the graph....")
df = CSV.read(
        HTTP.get("https://gist.githubusercontent.com/chriddyp/" *
        "cb5392c35661370d95f300086accea51/raw/" *
        "8e0768211f6b747c0db42a9ce9a0937dafcbd8b2/" *
        "indicators.csv").body
    )
df = dropmissing(df)
available_indicators = unique(df[:,Symbol("Indicator Name")])

app = Dash("Dash Layout", external_stylesheets=external_stylesheets) do
    html_div() do
        html_div() do
            html_div(style = (width="48%", display="inline-block")) do
                dcc_dropdown(
                    id = "xaxis-column",
                    options = [(label = i, value = i) for i in available_indicators],
                    value="Fertility rate, total (births per woman)"
                ),
                dcc_radioitems(
                    id="xaxis-type",
                    options=[(label = i, value = i)  for i in ["Linear", "Log"]],
                    value="Linear",
                    labelStyle=(display = "inline-block",)
                )
            end,
            html_div(style = (width="48%", display="inline-block", float="right")) do
                dcc_dropdown(
                    id = "yaxis-column",
                    options = [(label = i, value = i) for i in available_indicators],
                    value="Life expectancy at birth, total (years)"
                ),
                dcc_radioitems(
                    id="yaxis-type",
                    options=[(label = i, value = i)  for i in ["Linear", "Log"]],
                    value="Linear",
                    labelStyle=(display = "inline-block",)
                )
            end

        end,
        html_div(style=(width = "49%", display = "inline-block", padding = "0 20")) do
            dcc_graph(id = "crossfilter-indicator-scatter")
        end,
        html_div(style=(width = "49%", display = "inline-block")) do
            dcc_graph(id = "x-time-series"),
            dcc_graph(id = "y-time-series")
        end,
        html_div(style=(width = "49%", padding = "0px 20px 20px 20px")) do
            dcc_slider(
                id = "year-slider",
                min = minimum(df[:,:Year]),
                max = maximum(df[:,:Year]),
                value = minimum(df[:,:Year]),
                marks = Dict([Symbol(v)=>Symbol(v) for v in unique(df[:,:Year])]),
                step = nothing
                )
        end
    end
end

callback!(app,
    callid"""xaxis-column.value, yaxis-column.value, xaxis-type.value, yaxis-type.value, year-slider.value
    => crossfilter-indicator-scatter.figure""") do xaxis_column_name, yaxis_column_name, xaxis_type, yaxis_type, year_value
    dff = df[df.Year .== year_value, :]
    Plot(
        dff[dff[Symbol("Indicator Name")] .== xaxis_column_name, :Value],
        dff[dff[Symbol("Indicator Name")] .== yaxis_column_name, :Value],
        Layout(
            xaxis_type = xaxis_type == "Linear" ? "linear" : "log",
            xaxis_title = xaxis_column_name,
            yaxis_title = yaxis_column_name,
            yaxis_type = yaxis_type == "Linear" ? "linear" : "log",
            hovermode = "closest",
            height=450
        ),
        kind = "scatter",
        text = dff[dff[Symbol("Indicator Name")] .== yaxis_column_name, Symbol("Country Name")],
        customdata = dff[dff[Symbol("Indicator Name")] .== yaxis_column_name, Symbol("Country Name")],
        mode = "markers",
        marker_size = 15,
        marker_opacity = 0.5,
        marker_line_width = 0.5,
        marker_line_color = "white"

    )

end

function create_time_series(dff, axis_type, title)
    Plot(
        dff[:, :Year],
        dff[:, :Value],
        Layout(
            yaxis_type = axis_type == "Linear" ? "linear" : "log",
            xaxis_showgrid = false,
            annotations = [
                attr(
                    x = 0,
                    y = 0.85,
                    xanchor = "left",
                    yanchor = "bottom",
                    xref = "paper",
                    yref = "paper",
                    showarrow = false,
                    align = "left",
                    bgcolor = "rgba(255, 255, 255, 0.5)",
                    text = title
                )
            ],
            height=225
        ),
        mode = "lines+markers"
    )
end

callback!(app, callid"""
                crossfilter-indicator-scatter.hoverData,
                xaxis-column.value,
                xaxis-type.value => x-time-series.figure
                """) do hover_data, xaxis_column_name, axis_type

    country_name = isnothing(hover_data) ? "" : hover_data.points[1].customdata
    dff = df[df[:,Symbol("Country Name")] .==  country_name, :]
    dff = dff[dff[:,Symbol("Indicator Name")] .==  xaxis_column_name, :]
    title = "$(country_name)
$(xaxis_column_name)"
    create_time_series(dff, axis_type, title)
end

callback!(app, callid"""
                crossfilter-indicator-scatter.hoverData,
                yaxis-column.value,
                yaxis-type.value => y-time-series.figure
                """) do hover_data, yaxis_column_name, axis_type

    country_name = isnothing(hover_data) ? "" : hover_data.points[1].customdata
    dff = df[df[:,Symbol("Country Name")] .==  country_name, :]
    dff = dff[dff[:,Symbol("Indicator Name")] .==  yaxis_column_name, :]

    create_time_series(dff, axis_type, yaxis_column_name)
end

handler = make_handler(app, debug = true)
println("started at localhost:$(parse(Int, ARGS[1]))")
HTTP.serve(handler,"0.0.0.0",parse(Int,ARGS[1]))


text:
using IterativeSolvers
using FactCheck

n = 10
m = 6
srand(1234321)

include("common.jl")

##################
# Linear solvers #
##################

facts("Stationary solvers") do
    for T in (Float32, Float64, Complex64, Complex128)
        context("$T") do

        A=convert(Matrix{T}, randn(n,n))
        T<:Complex && (A+=convert(Matrix{T}, im*randn(n,n)))
        A+=A'
        A+=n*maximum(abs(A[:]))*diagm(ones(T,n)) #Diagonally dominant
        b=convert(Vector{T}, randn(n))
        T<:Complex && (b+=convert(Vector{T}, im*randn(n)))
        b=b/norm(b)
        x0=convert(Vector{T}, randn(n))
        T<:Complex && (x0+=convert(Vector{T}, im*randn(n)))
        x = A\b
        for solver in [jacobi, gauss_seidel]
            xi, ci=solver(A, b, maxiter=n^4)
            @fact ci.isconverged --> true
            @fact norm(x-xi) --> less_than(n^3*eps(typeof(real(b[1]))))
        end
        for solver in [jacobi!, gauss_seidel!]
            xi, ci=solver(copy(x0), A, b, maxiter=n^4)
            @fact ci.isconverged --> true
            @fact norm(x-xi) --> less_than(n^3*eps(typeof(real(b[1]))))
        end

        ω = 0.5
        for solver in [sor, ssor]
            xi, ci=solver(A, b, ω, maxiter=n^4)
            @fact ci.isconverged --> true
            @fact norm(x-xi) --> less_than(n^3*eps(typeof(real(b[1]))))
        end
        for solver in [sor!, ssor!]
            xi, ci=solver(copy(x0), A, b, ω, maxiter=n^4)
            @fact ci.isconverged --> true
            @fact norm(x-xi) --> less_than(n^3*eps(typeof(real(b[1]))))
        end

        end
    end
end

#Conjugate gradients
include("cg.jl")

#GMRES
facts("gmres") do

for T in (Float32, Float64, Complex64, Complex128)
    context("Matrix{$T}") do

    A = convert(Matrix{T}, randn(n,n))
    L = convert(Matrix{T}, randn(n,n))
    R = convert(Matrix{T}, randn(n,n))
    b = convert(Vector{T}, randn(n))
    if T <: Complex
        A += im*convert(Matrix{T}, randn(n,n))
        L += im*convert(Matrix{T}, randn(n,n))
        R += im*convert(Matrix{T}, randn(n,n))
        b += im*convert(Vector{T}, randn(n))
    end
    F = lufact(A)
    b = b/norm(b)

    x_gmres, c_gmres = gmres(A, b, L, R)
    @fact c_gmres.isconverged --> true
    @fact norm(A*x_gmres - b) --> less_than(√eps(real(one(T))))

    x_gmres, c_gmres = gmres(A, b, F; maxiter=1, restart=1)
    @fact c_gmres.isconverged --> true
    @fact norm(A*x_gmres - b) --> less_than(√eps(real(one(T))))

    x_gmres, c_gmres = gmres(A, b, 1, F; maxiter=1, restart=1)
    @fact c_gmres.isconverged --> true
    @fact norm(A*x_gmres - b) --> less_than(√eps(real(one(T))))
    end
end

for T in (Float64, Complex128)
    context("SparseMatrixCSC{$T}") do
    A = sprandn(n,n,0.5)+0.001*eye(T,n,n)
    L = sprandn(n,n,0.5)
    R = sprandn(n,n,0.5)
    b = convert(Vector{T}, randn(n))
    if T <: Complex
        A += im*sprandn(n,n,0.5)
        L += im*sprandn(n,n,0.5)
        R += im*sprandn(n,n,0.5)
        b += im*randn(n)
    end
    F = lufact(A)
    b = b / norm(b)

    x_gmres, c_gmres= gmres(A, b, L, R)
    @fact c_gmres.isconverged --> true
    @fact norm(A*x_gmres - b) --> less_than(√eps(real(one(T))))

    x_gmres, c_gmres = gmres(A, b, F; maxiter=1, restart=1)
    @fact c_gmres.isconverged --> true
    @fact norm(A*x_gmres - b) --> less_than(√eps(real(one(T))))

    x_gmres, c_gmres = gmres(A, b, 1, F; maxiter=1, restart=1)
    @fact c_gmres.isconverged --> true
    @fact norm(A*x_gmres - b) --> less_than(√eps(real(one(T))))
    end
end
end

#IDRs
facts("idrs") do

for T in (Float32, Float64, Complex64, Complex128)
    # without smoothing
    context("Matrix{$T}, without residual smoothing") do

    A = convert(Matrix{T}, randn(n,n))
    b = convert(Vector{T}, randn(n))
    if T <: Complex
        A += im*convert(Matrix{T}, randn(n,n))
        b += im*convert(Vector{T}, randn(n))
    end
    b = b/norm(b)

    x_idrs, c_idrs = idrs(A, b)
    @fact c_idrs.isconverged --> true
    @fact norm(A*x_idrs - b) --> less_than(√eps(real(one(T))))
    end

    # with smoothing
    context("Matrix{$T}, with residual smoothing") do

    A = convert(Matrix{T}, randn(n,n))
    b = convert(Vector{T}, randn(n))
    if T <: Complex
        A += im*convert(Matrix{T}, randn(n,n))
        b += im*convert(Vector{T}, randn(n))
    end
    b = b/norm(b)

    x_idrs, c_idrs = idrs(A, b; smoothing=true)
    @fact c_idrs.isconverged --> true
    @fact norm(A*x_idrs - b) --> less_than(√eps(real(one(T))))
    end
end

for T in (Float64, Complex128)
    context("SparseMatrixCSC{$T}") do
    A = sprandn(n,n,0.5)+0.001*eye(T,n,n)
    b = convert(Vector{T}, randn(n))
    if T <: Complex
        A += im*sprandn(n,n,0.5)
        b += im*randn(n)
    end
    b = b / norm(b)

    x_idrs, c_idrs= idrs(A, b)
    @fact c_idrs.isconverged --> true
    @fact norm(A*x_idrs - b) --> less_than(√eps(real(one(T))))
    end
end
end

facts("Chebyshev") do
for T in (Float32, Float64, Complex64, Complex128)
    context("Matrix{$T}") do
    A=convert(Matrix{T}, randn(n,n))
    T<:Complex && (A+=convert(Matrix{T}, im*randn(n,n)))
    A=A+A'
    A=A'*A #Construct SPD matrix
    b=convert(Vector{T}, randn(n))
    T<:Complex && (b+=convert(Vector{T}, im*randn(n)))
    b=b/norm(b)
    tol = 0.1 #For some reason Chebyshev is very slow
    v = eigvals(A)
    mxv = maximum(v)
    mnv = minimum(v)
    x_cheby, c_cheby= chebyshev(A, b, mxv+(mxv-mnv)/100, mnv-(mxv-mnv)/100, tol=tol, maxiter=10^5)
    @fact c_cheby.isconverged --> true
    @fact norm(A*x_cheby-b) --> less_than(tol)
    end
end
end

#######################
# Eigensystem solvers #
#######################

facts("simple eigensolvers") do
for T in (Float32, Float64, Complex64, Complex128)
    context("Matrix{$T}") do
    A=convert(Matrix{T}, randn(n,n))
    T<:Complex && (A+=convert(Matrix{T}, im*randn(n,n)))
    A=A+A' #Symmetric/Hermitian

    tol = (eltype(T) <: Complex ?2:1)*n^2*cond(A)*eps(real(one(T)))
    v = eigvals(A)

    ## Simple methods

    context("Power iteration") do
    eval_big = maximum(v) > abs(minimum(v)) ? maximum(v) : minimum(v)
    eval_pow = eigvals_power(A; tol=sqrt(eps(real(one(T)))), maxiter=2000)[1].val
    @fact norm(eval_big-eval_pow) --> less_than(tol)
    end

    context("Inverse iteration") do
    irnd = ceil(Int, rand()*(n-2))
    eval_rand = v[1+irnd] #Pick random eigenvalue
    # Perturb the eigenvalue by < 1/4 of the distance to the nearest eigenvalue
    eval_diff = min(abs(v[irnd]-eval_rand), abs(v[irnd+2]-eval_rand))
    σ = eval_rand + eval_diff/2*(rand()-.5)
    eval_ii = eigvals_ii(A, σ; tol=sqrt(eps(real(one(T)))), maxiter=2000)[1].val
    @fact norm(eval_rand-eval_ii) --> less_than(tol)
    end

    #context("Rayleigh quotient iteration") do
    #XXX broken?
    #l = eigvals_rqi(A, eigvals_rand, 2000, √eps()).val
    #@fact norm(eigvals_rand-l) --> less_than(tol)
    #end

    end
end
end

#Lanczos methods

type MyOp{T}
    buf::Matrix{T}
end

import Base: *, size, eltype

*(A::MyOp, x::AbstractVector) = A.buf*x
size(A::MyOp, i) = size(A.buf, i)
eltype(A::MyOp) = eltype(A.buf)

facts("eigvals_lanczos") do
for T in (Float32, Float64)
    context("Matrix{$T}") do
    A = convert(Matrix{T}, randn(n,n))
    A = A + A' #Symmetric
    v = eigvals(A)

    eval_lanczos, c_lanczos = eigvals_lanczos(A)
    @fact c_lanczos.isconverged --> true
    @fact norm(v - eval_lanczos) --> less_than(√eps(T))
    end

    context("Op{$T}") do
        A = MyOp(convert(Matrix{T}, randn(5,5)) |> t -> t + t')
        v = eigvals(Symmetric(A.buf))
        eval_lanczos, c_lanczos = eigvals_lanczos(A)
        @fact c_lanczos.isconverged --> true
        @fact norm(v - eval_lanczos) --> less_than(√eps(T))
    end
end
end

#Golub-Kahan-Lanczos singular values computation
include("svdl.jl")

include("lsqr.jl")
include("lsmr.jl")

#Randomized algorithms
include("rlinalg.jl")
include("rsvd.jl")
include("rsvd_fnkz.jl")

#Expensive tests - don't run by default
#include("matrixmarket.jl")
#include("matrixcollection.jl")

exitstatus()


text:
export SGD, Nesterov

export LearningRatePolicy, MomentumPolicy, SGDSolverParameters

@compat abstract type LearningRatePolicy end
@compat abstract type MomentumPolicy end


text:
__precompile__(true)
module Dolang

using DataStructures
using StaticArrays

const ARITH_SYMBOLS = Set([:+, :-, :*, :/, :^])
const DOLANG_FUNCTIONS = Set([:sin, :cos, :tan, :exp, :log, :log10])
push!(DOLANG_FUNCTIONS, ARITH_SYMBOLS...)

import Base: ==
import REPL

using Printf
export make_function, Der, FlatFunctionFactory, stringify, is_time_shift, time_shift,
       subs, csubs, steady_state, list_symbols, list_symbols!, list_variables,
       list_parameters

# come convenience methods
_replace_star_star(s::AbstractString) = replace(s, "**" => "^")

_to_expr(x::Expr) = x
_to_expr(x::Union{Symbol,Number}) = Expr(:block, x)
_to_expr(x::AbstractString) = _to_expr(Meta.parse(_replace_star_star(x)))

import SymEngine
deriv(eq::SymEngine.Basic, x) = SymEngine.diff(eq, x)
# somebody will probably object to it:
deriv(eq::Expr, x::Symbol) = Meta.parse(string(SymEngine.diff(SymEngine.Basic(eq), SymEngine.Basic(x))))
deriv(a::Symbol, b::Symbol) = a==b ? 1.0 : 0.0
@inline prep_deriv(eq) = SymEngine.Basic(eq)
@inline post_deriv(eq) = SymEngine.walk_expression(eq)

abstract type COrder end
abstract type FOrder end

# core files
include("symbolic.jl")
include("incidence.jl")
include("factory.jl")
include("util.jl")
include("compiler.jl")
include("printing.jl")
include("minilanguage.jl")
include("grammar.jl")

end  # module

max_issues_repo_name: EconForge/Dolang


text:
using Test
using WaveProp
using WaveProp.Nystrom
using WaveProp.Geometry
using WaveProp.Integration
using WaveProp.Mesh

@testset "Potential test" begin
    pde  = Helmholtz(;dim=3,k=1)
    Ω,M  = WaveProp.IO.gmsh_sphere(dim=2)
    Γ    = boundary(Ω)
    mesh = NystromMesh(view(M,Γ),order=1)
    𝓢    = SingleLayerPotential(pde,mesh)
    𝓓    = DoubleLayerPotential(pde,mesh)
    @test Nystrom.kernel_type(𝓢) == Nystrom.SingleLayer()
    @test Nystrom.kernel_type(𝓓) == Nystrom.DoubleLayer()
    σ    = γ₀(x->0.,mesh)
    u(x) = 𝓢[σ](x)
    x₀ = SVector(1.,1.,1.)
    @test u(x₀) ≈ 0
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp