text:
@testset "interface test" begin
    @testset "control gate" begin
        @test control(X) == ConstGate(control(LinQu.X_kernel))
        @test control(Rx)[π] == control(Rx[π])
    end
end


text:
###############################################################################
#
# An example of how to create autofilters with XlsxWriter.
#
# An autofilter is a way of adding drop down lists to the headers of a 2D
# range of worksheet data. This allows users to filter the data based on
# simple criteria so that some data is shown and some is hidden.
#
# Original Python Copyright 2013-2016, John McNamara, jmcnamara@cpan.org
# https://github.com/jmcnamara/XlsxWriter
#
#

using XlsxWriter
using DelimitedFiles

function test()

	wb = Workbook("autofilter.xlsx")

	# Add a worksheet for each autofilter example.
	sheet = Worksheet[]
	for i in 1:6
		push!(sheet, add_worksheet!(wb))
	end

	# Add a bold format for the headers.
	bold = add_format!(wb, Dict("bold"=>true))

	# Open a text file with autofilter example data.
	(data, headers) = readdlm("autofilter_data.txt", header=true)

	# Set up several sheets with the same data.
	for ws in sheet
		# Make the columns wider.
		set_column!(ws, "A:D", 12)
		# Make the header row larger and bold.
		set_row!(ws, 0, height=20, fmt=bold)
		# Add header text
		write_row!(ws, 0, 0, headers)
	end

	###############################################################################
	# Example 1. Autofilter without conditions.

	# Set the autofilter.
	autofilter!(sheet[1], "A1:D51")
	write_matrix!(sheet[1], 1, 0, data)

	###############################################################################
	# Example 2. Autofilter with a filter condition in the first column.

	# Autofilter range using Row-Column notation.
	autofilter!(sheet[2], 0, 0, 50, 3)

	# Add filter criteria. The placeholder "Region" in the filter is
	# ignored and can be any string that adds clarity to the expression.
	filter_column!(sheet[2], 0, "Region == East")

	# Hide the rows that don't match the filter criteria.

	for r in 1:size(data, 1)
		# Check for rows that match the filter.
		if data[r,1] != "East"
			set_row!(sheet[2], r, options=Dict("hidden"=>"1"))
		end
		write_row!(sheet[2], r, 0, data[r,:])
	end

	###############################################################################
	# Example 3. Autofilter with a dual filter condition in one of the columns.

	# Set the autofilter.
	autofilter!(sheet[3], "A1:D51")

	# Add filter criteria.
	filter_column!(sheet[3], "A", "x == East or x == South")

	# Hide the rows that don't match the filter criteria.
	for r in 1:size(data, 1)
		# Check for rows that match the filter.
		if data[r,1] != "East" && data[r,1] != "South"
			set_row!(sheet[3], r, options=Dict("hidden"=>true))
		end
		write_row!(sheet[3], r, 0, data[r,:])
	end

	###############################################################################
	# Example 4. Autofilter with filter conditions in two columns.

	# Set the autofilter.
	autofilter!(sheet[4], "A1:D51")

	# Add filter criteria.
	filter_column!(sheet[4], "A", "x == East")
	filter_column!(sheet[4], "C", "x > 3000 and x < 8000")

	# Hide the rows that don't match the filter criteria.
	for r in 1:size(data, 1)
		# Hide rows that don't match the filter.
		if !(data[r,1] == "East" && data[r,3] > 3000 && data[r,3] < 8000)
			set_row!(sheet[4], r, options=Dict("hidden"=>true))
		end
		write_row!(sheet[4], r, 0, data[r,:])
	end

	###############################################################################
	# Example 5. Autofilter with filter for blanks.
	# Create a blank cell in our test data.

	# Set the autofilter.
	autofilter!(sheet[5], "A1:D51")

	# Add filter criteria.
	filter_column!(sheet[5], "A", "x == Blanks")

	# Simulate a blank cell in the data.
	data[5, 1] = ""

	for r in 1:size(data, 1)
		# Check for rows that match the filter.
		if data[r,1] != ""
			# Row matches the filter, no further action required.
			set_row!(sheet[5], r, options=Dict("hidden"=>true))
		end
		write_row!(sheet[5], r, 0, data[r,:])
	end

	###############################################################################
	# Example 6. Autofilter with filter for non-blanks.
	# Set the autofilter.
	autofilter!(sheet[6], "A1:D51")

	# Add filter criteria.
	filter_column!(sheet[6], "A", "x == NonBlanks")

	for r in 1:size(data, 1)
		# Check for rows that match the filter.
		if data[r,1] == ""
			# Row matches the filter, no further action required.
			set_row!(sheet[6], r, options=Dict("hidden"=>true))
		end
		write_row!(sheet[6], r, 0, data[r,:])
	end

	close(wb)
	isfile("autofilter.xlsx")
end

test()

ext: jl


text:
using HorizonSideRobots
include("roblib.jl")

function main(r)
    num_steps_max=1
    side=Nord
    while ismarker(r)==false
        for _ in 1:2
            find_marker(r,side,num_steps_max)
            side=next(side)
        end
        num_steps_max+=1
    end
end

function find_marker(r,side,num_steps_max)
    for _ in 1:num_steps_max
        if ismarker(r)
            return nothing
        end
        move!(r,side)
    end
end

next(side::HorizonSide)=HorizonSide(mod(Int(side)+1,4))

r=Robot("9.sit", animate=true)

max_stars_repo_name: Savko-fokin/mirea-progs

max_issues_repo_name: Savko-fokin/mirea-progs

max_forks_repo_name: Savko-fokin/mirea-progs


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

mutable struct ConnectivityHop2 <: SwaggerModel
    type::Any # spec type: Union{ Nothing, String } # spec name: type
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    address::Any # spec type: Union{ Nothing, String } # spec name: address
    resourceId::Any # spec type: Union{ Nothing, String } # spec name: resourceId
    nextHopIds::Any # spec type: Union{ Nothing, Vector{String} } # spec name: nextHopIds
    previousHopIds::Any # spec type: Union{ Nothing, Vector{String} } # spec name: previousHopIds
    links::Any # spec type: Union{ Nothing, Vector{HopLink2} } # spec name: links
    previousLinks::Any # spec type: Union{ Nothing, Vector{HopLink2} } # spec name: previousLinks
    issues::Any # spec type: Union{ Nothing, Vector{ConnectivityIssue2} } # spec name: issues

    function ConnectivityHop2(;type=nothing, id=nothing, address=nothing, resourceId=nothing, nextHopIds=nothing, previousHopIds=nothing, links=nothing, previousLinks=nothing, issues=nothing)
        o = new()
        validate_property(ConnectivityHop2, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        validate_property(ConnectivityHop2, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(ConnectivityHop2, Symbol("address"), address)
        setfield!(o, Symbol("address"), address)
        validate_property(ConnectivityHop2, Symbol("resourceId"), resourceId)
        setfield!(o, Symbol("resourceId"), resourceId)
        validate_property(ConnectivityHop2, Symbol("nextHopIds"), nextHopIds)
        setfield!(o, Symbol("nextHopIds"), nextHopIds)
        validate_property(ConnectivityHop2, Symbol("previousHopIds"), previousHopIds)
        setfield!(o, Symbol("previousHopIds"), previousHopIds)
        validate_property(ConnectivityHop2, Symbol("links"), links)
        setfield!(o, Symbol("links"), links)
        validate_property(ConnectivityHop2, Symbol("previousLinks"), previousLinks)
        setfield!(o, Symbol("previousLinks"), previousLinks)
        validate_property(ConnectivityHop2, Symbol("issues"), issues)
        setfield!(o, Symbol("issues"), issues)
        o
    end
end # type ConnectivityHop2

const _property_map_ConnectivityHop2 = Dict{Symbol,Symbol}(Symbol("type")=>Symbol("type"), Symbol("id")=>Symbol("id"), Symbol("address")=>Symbol("address"), Symbol("resourceId")=>Symbol("resourceId"), Symbol("nextHopIds")=>Symbol("nextHopIds"), Symbol("previousHopIds")=>Symbol("previousHopIds"), Symbol("links")=>Symbol("links"), Symbol("previousLinks")=>Symbol("previousLinks"), Symbol("issues")=>Symbol("issues"))
const _property_types_ConnectivityHop2 = Dict{Symbol,String}(Symbol("type")=>"String", Symbol("id")=>"String", Symbol("address")=>"String", Symbol("resourceId")=>"String", Symbol("nextHopIds")=>"Vector{String}", Symbol("previousHopIds")=>"Vector{String}", Symbol("links")=>"Vector{HopLink2}", Symbol("previousLinks")=>"Vector{HopLink2}", Symbol("issues")=>"Vector{ConnectivityIssue2}")
Base.propertynames(::Type{ ConnectivityHop2 }) = collect(keys(_property_map_ConnectivityHop2))
Swagger.property_type(::Type{ ConnectivityHop2 }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ConnectivityHop2[name]))}
Swagger.field_name(::Type{ ConnectivityHop2 }, property_name::Symbol) =  _property_map_ConnectivityHop2[property_name]

function check_required(o::ConnectivityHop2)
    true
end

function validate_property(::Type{ ConnectivityHop2 }, name::Symbol, val)
end


text:
function FFtests()
  printstyled(" Quadratic Fields\n", color = :red)
  @time begin
    lq = fields(2, 1, fmpz(10)^4)
    @assert length(lq) == 6086

    lqreal = fields(2, 1, fmpz(10)^4, only_real = true)
    @assert length(lqreal) == 3043

    lq2 = [x for x in lq if istotally_real(x.field)]
    @assert length(lq2) == 3043
  end

  printstyled("\n Biquadratic Fields\n", color = :red)
  @time begin
    l = fields(4, 2, fmpz(10)^6)
    @assert length(l) ==  1014

    l1 = fields(4, 2, fmpz(10)^6, only_real = true)
    @assert length(l1) == 196

    l2 = [x for x in l if istotally_real(x.field)]
    @assert length(l2) == 196
  end

  printstyled("\n Direct product decomposition\n", color = :red)
  @time begin
    l_direct_product = fields(12, 4, fmpz(10)^13)
    @assert length(l_direct_product) == 25

    l_without = fields(12, 4, fmpz(10)^13, using_direct_product = false)
    @assert length(l_without) == 25

    d1 = length(fields(24, 13, fmpz(10)^24, simplify = false))
    d2 = length(fields(24, 13, fmpz(10)^24, using_direct_product = false, simplify = false))
    @assert d1 == d2
  end

  printstyled("\n Obstructions: degree 2\n", color = :red)
  @time begin
    G = GAP.Globals.SmallGroup(8, 4)
    L = GAP.Globals.DerivedSeries(G)
    forQ8_old = Hecke.check_obstruction(l, L, 2, [2])
    @assert length(forQ8_old) == 53
    forQ8 = [FieldsTower(x.field, x.generators_of_automorphisms, x.subfields) for x in forQ8_old]

    lQ8 = fields(8, 4, fmpz(10)^12)
    @assert length(lQ8) == 2

    lQ8_2 = fields(8, 4, forQ8, fmpz(10)^12)
    @assert length(lQ8_2) == 2

    lQ8real = fields(8, 4, fmpz(10)^12, only_real = true)
    @assert length(lQ8real) == 1
  end

  printstyled("\n Obstructions: prime degree\n", color = :red)
  @time begin
    @assert length(fields(54, 5, fmpz(10)^70, simplify = false)) == 3
    @assert length(fields(27, 3, fmpz(10)^60, simplify = false)) == 10
  end

  printstyled("\n Obstructions: prime_power_case\n", color = :red)
  @time begin
    l = fields(4, 2, fmpz(10)^5);
    G = GAP.Globals.SmallGroup(16, 9)
    L = GAP.Globals.DerivedSeries(G)
    lsieved = check_obstruction(l, L, 2, [4])
    @assert length(lsieved) == 78
    l = fields(4, 2, fmpz(10)^6);
    G = GAP.Globals.SmallGroup(32, 20)
    L = GAP.Globals.DerivedSeries(G)
    lsieved = check_obstruction(l, L, 2, [8])
    @assert length(lsieved) == 232
  end

  printstyled("\n Obstructions: reduction to cyclic case\n", color = :red)
  @time begin
    l = fields(8, 2, fmpz(10)^10)
    G = GAP.Globals.SmallGroup(96, 13)
    L = GAP.Globals.DerivedSeries(G)
    lsieved = check_obstruction(l, L, 2, [2, 6])
    @assert length(l) == length(lsieved)
    @assert length(fields(32, 30, fmpz(10)^39, simplify = false)) == 2
  end

  printstyled("\n Some examples\n", color = :red)
  @time begin
    @assert length(fields(6, 1, fmpz(10)^8, simplify = false)) == 107
    @assert length(fields(24, 8, fmpz(10)^30, simplify = false)) == 15
  end

  printstyled("\n Tests were successful", color = :red)
end


text:
module Symbolic
using ModelingToolkit
import ModelingToolkit.Constant

struct Lag{T,N}
    variable::T
    lag::N
end

#you can't register a constructor so we have to introduce this intermediate.
lag(x,i) = Lag(x,i)

@register lag(x,i)

@variables t x y z

function apply(expr::Equation, data)
    rhs = expr.rhs
    apply(rhs, data)
end

function apply(expr::Constant, data)
    # constants don't have an op field they are just a value.
    return expr.value
end

function apply(expr::Operation, data)
    # this method only exists to harmonize the API for Equation, Constant, and Operation
    # all the real work is happening in the three argument version below.
    apply(expr.op, expr, data)
end

# this uses the operation function as a trait, so that we can dispatch on it;
# allowing client code to extend the language using Multiple Dispatch.
function apply(op::Function, expr::Operation, data)
    # handles the case where there are no more arguments to find.
    # we assume this is a leaf node in the expression, which refers to a field in the data
    if length(expr.args) == 0
        return getproperty(data, expr.op.name)
    end
    anses = map(expr.args) do a
        apply(a, data)
    end
    return op(anses...)
end

function apply(op::typeof(lag), expr::Operation, data)
    var = expr.args[1].op.name
    sft = expr.args[2].value
    return getproperty(data[sft], var)
end

using Test
@testset "Expression Construction" begin
    eqn1 = z ~ x+y
    eqn2 = z ~ x + lag(x, 1)
    wdw = y ~ lag(x, -1) + lag(x, 0) + lag(x, 1)
end

@testset "Expression Evaluation" begin
@test apply(z~x+y, (x=1, y=2)) == 3
@test apply(z~x+y^2, (x=1, y=2)) == 5
@test apply(z~2*x+y^2, (x=1, y=2)) == 6
end

eqn2 = z ~ 1 + lag(x, 0)
# @show apply(eqn2, (x=(1:5), y=(1:4)))

# how do we get the formula to execute like a time series model where given a
# DataFrame with a "time colum" or a Vector{NamedTuple} with a timestamp field
# or a Vector{(t=Time(...), data=(...))}? I think we need to pass a context
# around from the initial Model through the recursion tree. We might a well make
# it a read/write context so that people can overload it however they want.

data = [
    (t=1, x=1),
    (t=2, x=2),
    (t=3, x=3),
    (t=4, x=4),
    (t=5, x=5),
]

struct LagData{T}
    data::T
    current::Int
end

import Base: getindex, getproperty

getindex(ld::LagData, i::Integer) = ld.data[ld.current + i]
function Base.getproperty(ld::LagData, s::Symbol)
    if s == :data || s == :current
        return getfield(ld, s)
    end
    getfield(ld.data[ld.current], s)
end

function tsapply(eqn, data, lags=(-1:1))
    center = ceil(Int, length(lags)/2)
    map(1-first(lags):length(data)-last(lags)) do i
        trip = LagData([data[i+l] for l in lags], center)
        apply(eqn, trip)
    end
end

# The tsapply approach, does not use a context, because it is able to define
# collection types that behave appropriately for the lag function and the
# regular apply primitive. LagData creates a container that you can index with a
# time offset like ld[-1].field to get the shifted version of field. Or you can
# index it without a time offset to get the center time (Δt=0). Then the
# function tsapply creates an outer loop that packages an arbitrary collection
# (assumed to contain rows indexed by time) into sliding windows.

# client code can extend the windowing by defining a new version of tsapply with
# a different type for the lags. We use a range to denote a symmetric window.

@testset "Basic Lag Formulas" begin
@test tsapply(z ~ 1 + lag(x, -1), data) == [2,3,4]
@test tsapply(z ~ 1 + lag(x,  0), data) == [3,4,5]
@test tsapply(z ~ 1 + lag(x, +1), data) == [4,5,6]

@test tsapply(z ~ 2 * lag(x, -1), data) == 2 .* (1:3)
@test tsapply(z ~ 2 * lag(x,  0), data) == 2 .* (2:4)
@test tsapply(z ~ 2 * lag(x, +1), data) == 2 .* (3:5)

@test tsapply(z ~ lag(x, -1) + lag(x, -1), data) == 2 .* (1:3)
@test tsapply(z ~ lag(x,  0) + lag(x,  0), data) == 2 .* (2:4)
@test tsapply(z ~ lag(x, +1) + lag(x, +1), data) == 2 .* (3:5)

@test tsapply(z ~ x + lag(x, -1), data) == [3,5,7]
@test tsapply(z ~ x + lag(x,  0), data) == [4,6,8]
@test tsapply(z ~ x + lag(x, +1), data) == [5,7,9]
end
end


text:
using Gadfly, RDatasets

set_default_plot_size(15cm, 15cm)

d = dataset("mlmRev", "Chem97")
idx = findall((d.Score .> 4) .& (d.Age .> 0))
plot(d[idx[1:4:end], :],
     x=:GCSEScore,
     ygroup=:Gender,
     xgroup=:Score,
     y=:Age,
     Geom.subplot_grid(Geom.point, free_y_axis=true),
     Theme(point_size=0.1mm, highlight_width=0.0mm))


text:
#
# Processing the weather HDF5 files in Julia 1.0
#
# This script opens an HDF5 file and creates DataFrame objects for the three
# datasets it contains. It also makes a plot of day/night
# average temperature from the "daily" dataset.
#

using HDF5
using Base.Filesystem
using DataFrames
using CategoricalArrays
using Gadfly
using Dates

#+
# Change this variables to point to the HDF file you want to process
#-
filePath = relpath("../data")
fileName = "monthly-201808.hdf5"

#+
# Open the hdf5 file
#-
h5f = h5open(joinpath(filePath,fileName),"r")

#+
# Extract references to the data and labes for the rawdata dataset
#-
dr0  = h5f["rawdata/block0_values"] # Floats dataset
dr1  = h5f["rawdata/block1_values"] # integer64 (datetime) dataset
dr0n = h5f["rawdata/block0_items"]
dr1n = h5f["rawdata/block1_items"]

#+
# Create DataFrames for the block0 (floats) and block1 (integer64 containing
# times in nanosecods). Note the tick used to transpose the read matrices.
# Build an array with the field labels from both datasets.
#-
dfr0 = DataFrame(read(dr0)')
dfr1 = DataFrame(read(dr1)')
drn  = vcat(read(dr1n),read(dr0n))

#+
# Concatenate both DataFrames and assign the column names.
#-
dfr = hcat(dfr1,dfr0)
names!(dfr,Symbol.(drn))

#+
# Convert the interger64 values to proper times. Remember those
# times are UTC, but the Julia datatime is plain (has no TZ info)
#-
dfr[:time] = unix2datetime.(div.(dfr[:time],1000000000))

#+
# Extract references to the data and labes for the hourly dataset
#-
dh0 =  h5f["hourly/block0_values"] # Floats dataset
dh1  = h5f["hourly/block1_values"] # integer dataset (count)
dh2  = h5f["hourly/block2_values"] # integer64 dataset (time)
dh0n = h5f["hourly/block0_items"]
dh1n = h5f["hourly/block1_items"]
dh2n = h5f["hourly/block2_items"]

#+
# Create dataframes from the dataset contents and prepare an
# array with the field labels
#
dfh0 = DataFrame(read(dh0)')
dfh1 = DataFrame(read(dh1)')
dfh2 = DataFrame(read(dh2)')
dfn = vcat(read(dh0n), read(dh1n), read(dh2n))

#+
# Concatenate the datasets, assign column names and convert int64 to datetime
#-
dfh = hcat(dfh0,dfh1,dfh2)
names!(dfh,Symbol.(dfn))
dfh[:time] = unix2datetime.(div.(dfh[:time],1000000000))

#+
# Get references for the daily datasets (contents and labels)
#-
dd0 = h5f["daily/block0_values"] # Integer64 (datetime)
dd1 = h5f["daily/block1_values"] # Floats
dd2 = h5f["daily/block2_values"] # Integers (count and daylight indicator)
dd0n = h5f["daily/block0_items"]
dd1n = h5f["daily/block1_items"]
dd2n = h5f["daily/block2_items"]

#+
# Create DataFrames and label array
#-
dfd0 = DataFrame(read(dd0)')
dfd1 = DataFrame(read(dd1)')
dfd2 = DataFrame(read(dd2)')
ddn = vcat(read(dd0n), read(dd1n), read(dd2n))

#+
# Concatenate dataframes, add column labels and convert datetime
#-
dfd = hcat(dfd0,dfd1,dfd2)
names!(dfd,Symbol.(ddn))
dfd[:date] = unix2datetime.(div.(dfd[:date],1000000000))

#+
# Convert the integer daytime indicator to a categorical, descriptive
# variable.
#-
dfd[:daytime] = CategoricalArray(map(x-> if x==1 return "Day" else return "Night" end , dfd[:daytime]))

#+
# Unstack the data, extracting just the temperature, so we have the day
# and night temperatures in the same row. Discard days with incomplete
# data and add a day-night difference column.
#-
dfdu = dropmissing!(unstack(dfd,[:date],:daytime,:temperature))
dfdu[:Difference] = dfdu[:Day] - dfdu[:Night]

#+
# Stack again the dataframe so we can use Gadfly to plot it.
#-
dfdt = stack(dfdu,[:Day,:Night,:Difference], variable_name=:daytime, value_name=:temperature)
plot(dfdt,x=:date,y=:temperature,color=:daytime,Geom.line)

max_stars_repo_name: jguillaumes/MeteoClient

max_issues_repo_name: jguillaumes/MeteoClient

max_forks_repo_name: jguillaumes/MeteoClient


text:
export BERLenv, RandomEnv, setup!, fitness, new_gen!, run!

abstract type BERLenv end

"Random BERL environment: returns random fitness"
mutable struct RandomEnv <: BERLenv
    name::String
    memory::Array # Array of storable data
end

"Sets up the environment (load data etc)"
function RandomEnv(cfg::Dict)
    RandomEnv("Random", [])
end

"Returns the fitness of one Cambrian Individual in this environment"
function fitness(indiv::Cambrian.Individual, env::BERLenv)
    rand(length(indiv.fitness))
end

"Sets up the env for the next generation"
function new_gen!(env::BERLenv, cfg::Dict=nothing)
    push!(env.memory, length(env.memory))
end

## Run function

function run!(e::Cambrian.Evolution, env::BERLenv)
    for i in (e.gen+1):e.cfg["n_gen"]
        step!(e)
        new_gen!(env, e.cfg)
    end
end


text:
# treereduce method without v0 is in Dagger.  Why does this live here?
function treereduce(f, xs, v0)
    length(xs) == 0 && return v0
    length(xs) == 1 && return xs[1]
    l = length(xs)
    m = div(l, 2)
    f(treereduce(f, xs[1:m]), treereduce(f, xs[m+1:end]))
end

dvcat(x...; dims) = vcat(x...)

function subtable(nds::NDSparse, r)
    NDSparse(keys(nds)[r], values(nds)[r], presorted=true, copy=false)
end

function subtable(t::IndexedTable, r)
    t[r]
end

function extrema_range(x::AbstractArray{T}, r::UnitRange) where T
    return first(x), last(x)
    if !(1 <= first(r) && last(r) <= length(x))
        throw(BoundsError(x, r))
    end

    isempty(r) && return extrema(x[r])
    mn = x[first(r)]
    mx = x[first(r)]
    @inbounds @simd for i in r
        mn = min(x[i], mn)
        mx = max(x[i], mx)
    end
    mn, mx
end

# Data loading utilities

Base.@deprecate load_table(args...;kwargs...) loadndsparse(args...; distributed=false, kwargs...)

function prettify_filename(f)
    f = basename(f)
    if endswith(lowercase(f), ".csv")
        f = f[1:end-4]
    end
    return f
end

function _loadtable_serial(T, file::Union{IO, AbstractString, AbstractArray};
                      delim=',',
                      indexcols=[],
                      datacols=nothing,
                      filenamecol=nothing,
                      agg=nothing,
                      presorted=false,
                      copy=false,
                      csvread=TextParse.csvread,
                      kwargs...)

    #println("LOADING ", file)
    count = Int[]

    samecols = nothing
    if indexcols !== nothing
        if indexcols isa Union{Int, Symbol}
            indexcols = (indexcols,)
        end
        samecols = collect(Iterators.filter(x->isa(x, Union{Tuple, AbstractArray}),
                                            indexcols))
    end
    if datacols !== nothing
        if datacols isa Union{Int, Symbol}
            datacols = (datacols,)
        end
        append!(samecols, collect(Iterators.filter(x->isa(x, Union{Tuple, AbstractArray}),
                                                   datacols)))
    end

    if samecols !== nothing
        samecols = map(x->map(string, x), samecols)
    end

    if isa(file, AbstractArray)
        cols, header, count = csvread(file, delim;
                                      samecols=samecols,
                                      kwargs...)
    else
        cols, header = csvread(file, delim; kwargs...)
    end

    header = map(string, header)

    if filenamecol !== nothing
        # mimick a file name column
        if filenamecol isa Pair
            filenamecol, f = filenamecol
        else
            f = prettify_filename
        end
        if isa(file, AbstractArray)
            namecol = reduce(vcat, fill.(f.(file), count))
        else
            namecol = fill(f(file), length(cols[1]))
        end

        cols = (namecol, cols...)
        if !isempty(header)
            pushfirst!(header, string(filenamecol))
        end
    end

    if isempty(cols)
        error("File contains no columns!")
    end

    n = length(first(cols))
    implicitindex = false

    ## Construct Index
    _indexcols = collect(map(x->lookupbyheader(header, x), indexcols))

    if isempty(_indexcols)
        implicitindex = true
        index = Columns(([1:n;],))
    else
        indexcolnames = map(indexcols, _indexcols) do name, i
            if i === nothing
                error("Cannot index by unknown column $name")
            else
                isa(name, Int) ? canonical_name(header[name]) : canonical_name(name)
            end
        end

        indexvecs = cols[_indexcols]

        nullableidx = findall(x->eltype(x) <: Union{DataValue,Nullable} || Missing <: eltype(x), indexvecs)
        if !isempty(nullableidx)
            badcol_names = header[_indexcols[nullableidx]]
            @warn("Indexed columns may contain Nullables or NAs. Column(s) with nullables: $(join(badcol_names, ", ", " and ")). This will result in wrong sorting.")
        end

        index = Columns(Tuple(indexvecs); names=indexcolnames)
    end

    ## Construct Data
    if datacols === nothing
        _datacols = setdiff(1:length(cols), _indexcols)
        datacols = header[_datacols]
    else
        _datacols = map(x->lookupbyheader(header, x), datacols)
    end

    if isempty(_datacols)
        error("""You must specify at least one data column.
                 Either all columns in the file were indexed, or datacols was explicitly set to an empty array.""")
    end

    datacolnames = map(datacols, _datacols) do name, i
        if i === nothing
            if isa(name, Int)
                error("Unknown column numbered $name specified in datacols")
            else
                return canonical_name(name) # use provided name for missing column
            end
        else
            isa(name, Int) ? canonical_name(header[name]) : canonical_name(name)
        end
    end

    datavecs = map(_datacols) do i
        if i === nothing
            # DataValueArray{Union{}}(n) # missing column
            fill(missing, n)
        else
            cols[i]
        end
    end

    data = Columns(Tuple(datavecs); names=datacolnames)

    if T<:IndexedTable && implicitindex
        table(data, copy = copy), true
    else
        convert(T, index, data, copy = copy), implicitindex
    end
end

function lookupbyheader(header, key)
    if isa(key, Symbol)
        return lookupbyheader(header, string(key))
    elseif isa(key, String)
        return findfirst(x->x==key, header)
    elseif isa(key, Int)
        return 0 < key <= length(header) ? key : nothing
    elseif isa(key, Tuple) || isa(key, Vector)
        for k in key
            x = lookupbyheader(header, k)
            x != 0 && return x
        end
        return nothing
    end
end

canonical_name(n::Symbol) = n
canonical_name(n::String) = Symbol(replace(n, r"\s" => "_"))
canonical_name(n::Union{Tuple, Vector}) = canonical_name(first(n))

function _repeated(x, n)
    Iterators.repeated(x,n)
end

function approx_size(cs::Columns)
    sum(map(approx_size, astuple(columns(cs))))
end

function approx_size(t::NDSparse)
    approx_size(t.data) + approx_size(t.index)
end

function approx_size(pa::PooledArray)
    approx_size(pa.refs) + approx_size(pa.pool) * 2
end

function approx_size(t::IndexedTable)
    approx_size(rows(t))
end

function approx_size(x::StringArray)
    approx_size(x.buffer) + approx_size(x.offsets) + approx_size(x.lengths)
end

# The following is not inferable, this is OK because the only place we use
# this doesn't need it.

function _map_params(f, T, S)
    (f(_tuple_type_head(T), _tuple_type_head(S)), _map_params(f, _tuple_type_tail(T), _tuple_type_tail(S))...)
end

_map_params(f, T::Type{Tuple{}},S::Type{Tuple{}}) = ()

map_params(f, ::Type{T}, ::Type{S}) where {T,S} = f(T,S)
@inline _tuple_type_head(::Type{T}) where {T<:Tuple} = Base.tuple_type_head(T)
@inline _tuple_type_tail(::Type{T}) where {T<:Tuple} = Base.tuple_type_tail(T)

#function map_params{N}(f, T::Type{T} where T<:Tuple{Vararg{Any,N}}, S::Type{S} where S<: Tuple{Vararg{Any,N}})
Base.@pure function map_params(f, ::Type{T}, ::Type{S}) where {T<:Tuple,S<:Tuple}
    if fieldcount(T) != fieldcount(S)
        MethodError(map_params, (typeof(f), T,S))
    end
    Tuple{_map_params(f, T,S)...}
end

_tuple_type_head(T::Type{NT}) where {NT<: NamedTuple} = fieldtype(NT, 1)

Base.@pure function _tuple_type_tail(T::Type{NamedTuple{n,ts}}) where {n,ts}
    _tuple_type_tail(ts)
end

Base.@pure @generated function map_params(f, ::Type{T}, ::Type{S}) where {T<:NamedTuple,S<:NamedTuple}
    if fieldnames(T) != fieldnames(S)
        MethodError(map_params, (T,S))
    end
    :(NamedTuple{$(fieldnames(T)), Tuple{_map_params(f, T, S)...}})
end

function randomsample(n, r::AbstractRange)
    k = 0
    taken = Set{eltype(r)}()
    output = eltype(r)[]
    while k < n && k < length(r)
        x = rand(r)
        if !(x in taken)
            push!(taken, x)
            push!(output, x)
            k += 1
        end
    end
    return output
end

function tuplesetindex(x::Tuple{Vararg{Any,N}}, v, i) where N
    ntuple(Val(N)) do j
        i == j ? v : x[j]
    end
end

@inline tuplesetindex(x::NamedTuple, v, i::Symbol) = (; x..., i => v)

@inline function tuplesetindex(x::NamedTuple{N}, v, i::Int) where N
    tuplesetindex(x, v, getfield(N, i))
end

function tuplesetindex(x::Union{NamedTuple, Tuple}, v::Tuple, i::Tuple)
    reduce((t, j)->tuplesetindex(t, v[j], i[j]), 1:length(i), init=x)
end


text:
"""
    SelectionMethod

Abstract type defining a method for choosing collocation points from spline
knots.
"""
abstract type SelectionMethod end

"""
    AvgKnots <: SelectionMethod

Each collocation point is chosen as a sliding average over `k - 1` knots:

```math
x_i = \\frac{1}{k - 1} ∑_{j = 1}^{k - 1} t_{i + j}
```

The resulting collocation points are sometimes called Greville sites
(de Boor 2001).
"""
struct AvgKnots <: SelectionMethod end

struct GrevilleSiteIterator{Basis <: AbstractBSplineBasis}
    basis :: Basis
end

BSplines.basis(it::GrevilleSiteIterator) = it.basis

# These are needed for `collect(it)`
Base.length(it::GrevilleSiteIterator) = length(basis(it))
Base.eltype(it::GrevilleSiteIterator) = float(eltype(knots(basis(it))))

@inline function Base.iterate(it::GrevilleSiteIterator, state = nothing)
    B = basis(it)
    N = length(B)
    T = eltype(it)

    if state === nothing
        i = 0
        compensation = zero(T)
    else
        i, tsum_prev, compensation = state
        i == N && return nothing
    end

    # For recombined bases, skip points at the boundaries.
    cl, cr = num_constraints(B)  # left/right constraints
    ii = i + cl

    k = order(B)
    ts = knots(B)
    lims = boundaries(B)

    if state === nothing
        tsum = zero(eltype(ts))
        for j = 2:k
            tsum += @inbounds ts[ii + j]
        end
    else
        # Optimised window averaging: reuse result from previous iteration.
        #
        # We use a two-step Kahan summation to avoid roundoff errors.
        # https://en.wikipedia.org/wiki/Kahan_summation_algorithm
        #
        # The non-compensated summation would be written as:
        # tsum = tsum_prev + ts[ii + k] - ts[ii + 1]

        # 1. Add next knot
        y = @inbounds ts[ii + k] - compensation
        tsum = tsum_prev + y
        compensation = (tsum - tsum_prev) - y

        tsum_prev = tsum

        # 2. Remove previous knot
        y = @inbounds -ts[ii + 1] - compensation
        tsum = tsum_prev + y
        compensation = (tsum - tsum_prev) - y
    end

    x = T(tsum / (k - 1))

    # If this is the first (last) point and there are no BCs on the left
    # (right), we make sure that the point exactly matches the left (right)
    # boundary.
    if cl === 0 && i == 0
        x = T(first(lims))
    elseif cr === 0 && i == N - 1
        x = T(last(lims))
    end

    # Make sure that the point is inside the domain.
    # This may not be the case if end knots have multiplicity less than k.
    x = clamp(x, lims...)

    x, (i + 1, tsum, compensation)
end

"""
    collocation_points(B::AbstractBSplineBasis; method=Collocation.AvgKnots())

Define and return adapted collocation points for evaluation of splines.

The number of returned collocation points is equal to the number of functions in
the basis.

Note that if `B` is a [`RecombinedBSplineBasis`](@ref) (adapted for boundary value
problems), collocation points are not included at the boundaries, since the
boundary conditions are implicitly satisfied by the basis.

In principle, the choice of collocation points is not unique.
The selection method can be chosen via the `method` argument.
For now, only a single method is accepted:

- [`Collocation.AvgKnots()`](@ref)

See also [`collocation_points!`](@ref).
"""
function collocation_points(
        B::AbstractBSplineBasis; method::SelectionMethod=AvgKnots())
    x = similar(knots(B), length(B))
    collocation_points!(x, B, method=method)
end

"""
    collocation_points!(x::AbstractVector, B::AbstractBSplineBasis;
                        method::SelectionMethod=AvgKnots())

Fill vector with collocation points for evaluation of splines.

See [`collocation_points`](@ref) for details.
"""
function collocation_points!(
        x::AbstractVector, B::AbstractBSplineBasis;
        method::SelectionMethod=AvgKnots())
    N = length(B)
    if N != length(x)
        throw(ArgumentError(
            "number of collocation points must match number of B-splines"))
    end
    _collocation_points!(x, B, method)
end

_collocation_points!(xs, B::AbstractBSplineBasis, ::AvgKnots) =
    copyto!(xs, GrevilleSiteIterator(B))

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
include("euler/euler.jl")
using .Primes: prime_numbers
using .Numbers: sum_proper_divisors
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(n::Int)::Int
    sum_factors = repeat([0], n)
    primes_numbers = prime_numbers(trunc(Int, √n))
    result = 0
    for i ∈ 1:n sum_factors[i] = trunc(Int, sum_proper_divisors(i, primes_numbers)) end
    for i ∈ 2:n
        j = sum_factors[i]
        if j ≠ i && j ≤ n && sum_factors[j] == i result += i end
    end
    return result
end

compute(1_000)

compute(10_000)

@benchmark compute(10_000)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sunlinsol/sunlinsol_band.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __SUNBandLinearSolver(y::N_Vector, A::SUNMatrix)
    ccall((:SUNBandLinearSolver, libsundials_sunlinsolband), SUNLinearSolver, (N_Vector, SUNMatrix), y, A)
end

function SUNBandLinearSolver(y, A)
    __y = convert(NVector, y)
    __SUNBandLinearSolver(convert(N_Vector, __y), A)
end

function SUNLinSolGetType_Band(S::SUNLinearSolver)
    ccall((:SUNLinSolGetType_Band, libsundials_sunlinsolband), SUNLinearSolver_Type, (SUNLinearSolver,), S)
end

function SUNLinSolInitialize_Band(S::SUNLinearSolver)
    ccall((:SUNLinSolInitialize_Band, libsundials_sunlinsolband), Cint, (SUNLinearSolver,), S)
end

function SUNLinSolSetup_Band(S::SUNLinearSolver, A::SUNMatrix)
    ccall((:SUNLinSolSetup_Band, libsundials_sunlinsolband), Cint, (SUNLinearSolver, SUNMatrix), S, A)
end

function __SUNLinSolSolve_Band(S::SUNLinearSolver, A::SUNMatrix, x::N_Vector, b::N_Vector, tol::realtype)
    ccall((:SUNLinSolSolve_Band, libsundials_sunlinsolband), Cint, (SUNLinearSolver, SUNMatrix, N_Vector, N_Vector, realtype), S, A, x, b, tol)
end

function SUNLinSolSolve_Band(S, A, x, b, tol)
    __x = convert(NVector, x)
    __b = convert(NVector, b)
    __SUNLinSolSolve_Band(S, A, convert(N_Vector, __x), convert(N_Vector, __b), tol)
end

function SUNLinSolLastFlag_Band(S::SUNLinearSolver)
    ccall((:SUNLinSolLastFlag_Band, libsundials_sunlinsolband), Clong, (SUNLinearSolver,), S)
end

function SUNLinSolSpace_Band(S::SUNLinearSolver, lenrwLS, leniwLS)
    ccall((:SUNLinSolSpace_Band, libsundials_sunlinsolband), Cint, (SUNLinearSolver, Ptr{Clong}, Ptr{Clong}), S, lenrwLS, leniwLS)
end

function SUNLinSolFree_Band(S::SUNLinearSolver)
    ccall((:SUNLinSolFree_Band, libsundials_sunlinsolband), Cint, (SUNLinearSolver,), S)
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sunlinsol/sunlinsol_dense.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __SUNDenseLinearSolver(y::N_Vector, A::SUNMatrix)
    ccall((:SUNDenseLinearSolver, libsundials_sunlinsoldense), SUNLinearSolver, (N_Vector, SUNMatrix), y, A)
end

function SUNDenseLinearSolver(y, A)
    __y = convert(NVector, y)
    __SUNDenseLinearSolver(convert(N_Vector, __y), A)
end

function SUNLinSolGetType_Dense(S::SUNLinearSolver)
    ccall((:SUNLinSolGetType_Dense, libsundials_sunlinsoldense), SUNLinearSolver_Type, (SUNLinearSolver,), S)
end

function SUNLinSolInitialize_Dense(S::SUNLinearSolver)
    ccall((:SUNLinSolInitialize_Dense, libsundials_sunlinsoldense), Cint, (SUNLinearSolver,), S)
end

function SUNLinSolSetup_Dense(S::SUNLinearSolver, A::SUNMatrix)
    ccall((:SUNLinSolSetup_Dense, libsundials_sunlinsoldense), Cint, (SUNLinearSolver, SUNMatrix), S, A)
end

function __SUNLinSolSolve_Dense(S::SUNLinearSolver, A::SUNMatrix, x::N_Vector, b::N_Vector, tol::realtype)
    ccall((:SUNLinSolSolve_Dense, libsundials_sunlinsoldense), Cint, (SUNLinearSolver, SUNMatrix, N_Vector, N_Vector, realtype), S, A, x, b, tol)
end

function SUNLinSolSolve_Dense(S, A, x, b, tol)
    __x = convert(NVector, x)
    __b = convert(NVector, b)
    __SUNLinSolSolve_Dense(S, A, convert(N_Vector, __x), convert(N_Vector, __b), tol)
end

function SUNLinSolLastFlag_Dense(S::SUNLinearSolver)
    ccall((:SUNLinSolLastFlag_Dense, libsundials_sunlinsoldense), Clong, (SUNLinearSolver,), S)
end

function SUNLinSolSpace_Dense(S::SUNLinearSolver, lenrwLS, leniwLS)
    ccall((:SUNLinSolSpace_Dense, libsundials_sunlinsoldense), Cint, (SUNLinearSolver, Ptr{Clong}, Ptr{Clong}), S, lenrwLS, leniwLS)
end

function SUNLinSolFree_Dense(S::SUNLinearSolver)
    ccall((:SUNLinSolFree_Dense, libsundials_sunlinsoldense), Cint, (SUNLinearSolver,), S)
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sunlinsol/sunlinsol_klu.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __SUNKLU(y::N_Vector, A::SUNMatrix)
    ccall((:SUNKLU, libsundials_sunlinsolklu), SUNLinearSolver, (N_Vector, SUNMatrix), y, A)
end

function SUNKLU(y, A)
    __y = convert(NVector, y)
    __SUNKLU(convert(N_Vector, __y), A)
end

function __SUNKLUReInit(S::SUNLinearSolver, A::SUNMatrix, nnz::sunindextype, reinit_type::Cint)
    ccall((:SUNKLUReInit, libsundials_sunlinsolklu), Cint, (SUNLinearSolver, SUNMatrix, sunindextype, Cint), S, A, nnz, reinit_type)
end

function SUNKLUReInit(S, A, nnz, reinit_type)
    __SUNKLUReInit(S, A, nnz, convert(Cint, reinit_type))
end

function __SUNKLUSetOrdering(S::SUNLinearSolver, ordering_choice::Cint)
    ccall((:SUNKLUSetOrdering, libsundials_sunlinsolklu), Cint, (SUNLinearSolver, Cint), S, ordering_choice)
end

function SUNKLUSetOrdering(S, ordering_choice)
    __SUNKLUSetOrdering(S, convert(Cint, ordering_choice))
end

function SUNLinSolGetType_KLU(S::SUNLinearSolver)
    ccall((:SUNLinSolGetType_KLU, libsundials_sunlinsolklu), SUNLinearSolver_Type, (SUNLinearSolver,), S)
end

function SUNLinSolInitialize_KLU(S::SUNLinearSolver)
    ccall((:SUNLinSolInitialize_KLU, libsundials_sunlinsolklu), Cint, (SUNLinearSolver,), S)
end

function SUNLinSolSetup_KLU(S::SUNLinearSolver, A::SUNMatrix)
    ccall((:SUNLinSolSetup_KLU, libsundials_sunlinsolklu), Cint, (SUNLinearSolver, SUNMatrix), S, A)
end

function __SUNLinSolSolve_KLU(S::SUNLinearSolver, A::SUNMatrix, x::N_Vector, b::N_Vector, tol::realtype)
    ccall((:SUNLinSolSolve_KLU, libsundials_sunlinsolklu), Cint, (SUNLinearSolver, SUNMatrix, N_Vector, N_Vector, realtype), S, A, x, b, tol)
end

function SUNLinSolSolve_KLU(S, A, x, b, tol)
    __x = convert(NVector, x)
    __b = convert(NVector, b)
    __SUNLinSolSolve_KLU(S, A, convert(N_Vector, __x), convert(N_Vector, __b), tol)
end

function SUNLinSolLastFlag_KLU(S::SUNLinearSolver)
    ccall((:SUNLinSolLastFlag_KLU, libsundials_sunlinsolklu), Clong, (SUNLinearSolver,), S)
end

function SUNLinSolSpace_KLU(S::SUNLinearSolver, lenrwLS, leniwLS)
    ccall((:SUNLinSolSpace_KLU, libsundials_sunlinsolklu), Cint, (SUNLinearSolver, Ptr{Clong}, Ptr{Clong}), S, lenrwLS, leniwLS)
end

function SUNLinSolFree_KLU(S::SUNLinearSolver)
    ccall((:SUNLinSolFree_KLU, libsundials_sunlinsolklu), Cint, (SUNLinearSolver,), S)
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sunlinsol/sunlinsol_pcg.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __SUNPCG(y::N_Vector, pretype::Cint, maxl::Cint)
    ccall((:SUNPCG, libsundials_sunlinsolpcg), SUNLinearSolver, (N_Vector, Cint, Cint), y, pretype, maxl)
end

function SUNPCG(y, pretype, maxl)
    __y = convert(NVector, y)
    __SUNPCG(convert(N_Vector, __y), convert(Cint, pretype), convert(Cint, maxl))
end

function __SUNPCGSetPrecType(S::SUNLinearSolver, pretype::Cint)
    ccall((:SUNPCGSetPrecType, libsundials_sunlinsolpcg), Cint, (SUNLinearSolver, Cint), S, pretype)
end

function SUNPCGSetPrecType(S, pretype)
    __SUNPCGSetPrecType(S, convert(Cint, pretype))
end

function __SUNPCGSetMaxl(S::SUNLinearSolver, maxl::Cint)
    ccall((:SUNPCGSetMaxl, libsundials_sunlinsolpcg), Cint, (SUNLinearSolver, Cint), S, maxl)
end

function SUNPCGSetMaxl(S, maxl)
    __SUNPCGSetMaxl(S, convert(Cint, maxl))
end

function SUNLinSolGetType_PCG(S::SUNLinearSolver)
    ccall((:SUNLinSolGetType_PCG, libsundials_sunlinsolpcg), SUNLinearSolver_Type, (SUNLinearSolver,), S)
end

function SUNLinSolInitialize_PCG(S::SUNLinearSolver)
    ccall((:SUNLinSolInitialize_PCG, libsundials_sunlinsolpcg), Cint, (SUNLinearSolver,), S)
end

function SUNLinSolSetATimes_PCG(S::SUNLinearSolver, A_data, ATimes::ATimesFn)
    ccall((:SUNLinSolSetATimes_PCG, libsundials_sunlinsolpcg), Cint, (SUNLinearSolver, Ptr{Cvoid}, ATimesFn), S, A_data, ATimes)
end

function SUNLinSolSetPreconditioner_PCG(S::SUNLinearSolver, P_data, Pset::PSetupFn, Psol::PSolveFn)
    ccall((:SUNLinSolSetPreconditioner_PCG, libsundials_sunlinsolpcg), Cint, (SUNLinearSolver, Ptr{Cvoid}, PSetupFn, PSolveFn), S, P_data, Pset, Psol)
end

function __SUNLinSolSetScalingVectors_PCG(S::SUNLinearSolver, s::N_Vector, nul::N_Vector)
    ccall((:SUNLinSolSetScalingVectors_PCG, libsundials_sunlinsolpcg), Cint, (SUNLinearSolver, N_Vector, N_Vector), S, s, nul)
end

function SUNLinSolSetScalingVectors_PCG(S, s, nul)
    __s = convert(NVector, s)
    __nul = convert(NVector, nul)
    __SUNLinSolSetScalingVectors_PCG(S, convert(N_Vector, __s), convert(N_Vector, __nul))
end

function SUNLinSolSetup_PCG(S::SUNLinearSolver, nul::SUNMatrix)
    ccall((:SUNLinSolSetup_PCG, libsundials_sunlinsolpcg), Cint, (SUNLinearSolver, SUNMatrix), S, nul)
end

function __SUNLinSolSolve_PCG(S::SUNLinearSolver, nul::SUNMatrix, x::N_Vector, b::N_Vector, tol::realtype)
    ccall((:SUNLinSolSolve_PCG, libsundials_sunlinsolpcg), Cint, (SUNLinearSolver, SUNMatrix, N_Vector, N_Vector, realtype), S, nul, x, b, tol)
end

function SUNLinSolSolve_PCG(S, nul, x, b, tol)
    __x = convert(NVector, x)
    __b = convert(NVector, b)
    __SUNLinSolSolve_PCG(S, nul, convert(N_Vector, __x), convert(N_Vector, __b), tol)
end

function SUNLinSolNumIters_PCG(S::SUNLinearSolver)
    ccall((:SUNLinSolNumIters_PCG, libsundials_sunlinsolpcg), Cint, (SUNLinearSolver,), S)
end

function SUNLinSolResNorm_PCG(S::SUNLinearSolver)
    ccall((:SUNLinSolResNorm_PCG, libsundials_sunlinsolpcg), realtype, (SUNLinearSolver,), S)
end

function SUNLinSolResid_PCG(S::SUNLinearSolver)
    ccall((:SUNLinSolResid_PCG, libsundials_sunlinsolpcg), N_Vector, (SUNLinearSolver,), S)
end

function SUNLinSolLastFlag_PCG(S::SUNLinearSolver)
    ccall((:SUNLinSolLastFlag_PCG, libsundials_sunlinsolpcg), Clong, (SUNLinearSolver,), S)
end

function SUNLinSolSpace_PCG(S::SUNLinearSolver, lenrwLS, leniwLS)
    ccall((:SUNLinSolSpace_PCG, libsundials_sunlinsolpcg), Cint, (SUNLinearSolver, Ptr{Clong}, Ptr{Clong}), S, lenrwLS, leniwLS)
end

function SUNLinSolFree_PCG(S::SUNLinearSolver)
    ccall((:SUNLinSolFree_PCG, libsundials_sunlinsolpcg), Cint, (SUNLinearSolver,), S)
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sunlinsol/sunlinsol_spbcgs.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __SUNSPBCGS(y::N_Vector, pretype::Cint, maxl::Cint)
    ccall((:SUNSPBCGS, libsundials_sunlinsolspbcgs), SUNLinearSolver, (N_Vector, Cint, Cint), y, pretype, maxl)
end

function SUNSPBCGS(y, pretype, maxl)
    __y = convert(NVector, y)
    __SUNSPBCGS(convert(N_Vector, __y), convert(Cint, pretype), convert(Cint, maxl))
end

function __SUNSPBCGSSetPrecType(S::SUNLinearSolver, pretype::Cint)
    ccall((:SUNSPBCGSSetPrecType, libsundials_sunlinsolspbcgs), Cint, (SUNLinearSolver, Cint), S, pretype)
end

function SUNSPBCGSSetPrecType(S, pretype)
    __SUNSPBCGSSetPrecType(S, convert(Cint, pretype))
end

function __SUNSPBCGSSetMaxl(S::SUNLinearSolver, maxl::Cint)
    ccall((:SUNSPBCGSSetMaxl, libsundials_sunlinsolspbcgs), Cint, (SUNLinearSolver, Cint), S, maxl)
end

function SUNSPBCGSSetMaxl(S, maxl)
    __SUNSPBCGSSetMaxl(S, convert(Cint, maxl))
end

function SUNLinSolGetType_SPBCGS(S::SUNLinearSolver)
    ccall((:SUNLinSolGetType_SPBCGS, libsundials_sunlinsolspbcgs), SUNLinearSolver_Type, (SUNLinearSolver,), S)
end

function SUNLinSolInitialize_SPBCGS(S::SUNLinearSolver)
    ccall((:SUNLinSolInitialize_SPBCGS, libsundials_sunlinsolspbcgs), Cint, (SUNLinearSolver,), S)
end

function SUNLinSolSetATimes_SPBCGS(S::SUNLinearSolver, A_data, ATimes::ATimesFn)
    ccall((:SUNLinSolSetATimes_SPBCGS, libsundials_sunlinsolspbcgs), Cint, (SUNLinearSolver, Ptr{Cvoid}, ATimesFn), S, A_data, ATimes)
end

function SUNLinSolSetPreconditioner_SPBCGS(S::SUNLinearSolver, P_data, Pset::PSetupFn, Psol::PSolveFn)
    ccall((:SUNLinSolSetPreconditioner_SPBCGS, libsundials_sunlinsolspbcgs), Cint, (SUNLinearSolver, Ptr{Cvoid}, PSetupFn, PSolveFn), S, P_data, Pset, Psol)
end

function __SUNLinSolSetScalingVectors_SPBCGS(S::SUNLinearSolver, s1::N_Vector, s2::N_Vector)
    ccall((:SUNLinSolSetScalingVectors_SPBCGS, libsundials_sunlinsolspbcgs), Cint, (SUNLinearSolver, N_Vector, N_Vector), S, s1, s2)
end

function SUNLinSolSetScalingVectors_SPBCGS(S, s1, s2)
    __s1 = convert(NVector, s1)
    __s2 = convert(NVector, s2)
    __SUNLinSolSetScalingVectors_SPBCGS(S, convert(N_Vector, __s1), convert(N_Vector, __s2))
end

function SUNLinSolSetup_SPBCGS(S::SUNLinearSolver, A::SUNMatrix)
    ccall((:SUNLinSolSetup_SPBCGS, libsundials_sunlinsolspbcgs), Cint, (SUNLinearSolver, SUNMatrix), S, A)
end

function __SUNLinSolSolve_SPBCGS(S::SUNLinearSolver, A::SUNMatrix, x::N_Vector, b::N_Vector, tol::realtype)
    ccall((:SUNLinSolSolve_SPBCGS, libsundials_sunlinsolspbcgs), Cint, (SUNLinearSolver, SUNMatrix, N_Vector, N_Vector, realtype), S, A, x, b, tol)
end

function SUNLinSolSolve_SPBCGS(S, A, x, b, tol)
    __x = convert(NVector, x)
    __b = convert(NVector, b)
    __SUNLinSolSolve_SPBCGS(S, A, convert(N_Vector, __x), convert(N_Vector, __b), tol)
end

function SUNLinSolNumIters_SPBCGS(S::SUNLinearSolver)
    ccall((:SUNLinSolNumIters_SPBCGS, libsundials_sunlinsolspbcgs), Cint, (SUNLinearSolver,), S)
end

function SUNLinSolResNorm_SPBCGS(S::SUNLinearSolver)
    ccall((:SUNLinSolResNorm_SPBCGS, libsundials_sunlinsolspbcgs), realtype, (SUNLinearSolver,), S)
end

function SUNLinSolResid_SPBCGS(S::SUNLinearSolver)
    ccall((:SUNLinSolResid_SPBCGS, libsundials_sunlinsolspbcgs), N_Vector, (SUNLinearSolver,), S)
end

function SUNLinSolLastFlag_SPBCGS(S::SUNLinearSolver)
    ccall((:SUNLinSolLastFlag_SPBCGS, libsundials_sunlinsolspbcgs), Clong, (SUNLinearSolver,), S)
end

function SUNLinSolSpace_SPBCGS(S::SUNLinearSolver, lenrwLS, leniwLS)
    ccall((:SUNLinSolSpace_SPBCGS, libsundials_sunlinsolspbcgs), Cint, (SUNLinearSolver, Ptr{Clong}, Ptr{Clong}), S, lenrwLS, leniwLS)
end

function SUNLinSolFree_SPBCGS(S::SUNLinearSolver)
    ccall((:SUNLinSolFree_SPBCGS, libsundials_sunlinsolspbcgs), Cint, (SUNLinearSolver,), S)
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sunlinsol/sunlinsol_spfgmr.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __SUNSPFGMR(y::N_Vector, pretype::Cint, maxl::Cint)
    ccall((:SUNSPFGMR, libsundials_sunlinsolspfgmr), SUNLinearSolver, (N_Vector, Cint, Cint), y, pretype, maxl)
end

function SUNSPFGMR(y, pretype, maxl)
    __y = convert(NVector, y)
    __SUNSPFGMR(convert(N_Vector, __y), convert(Cint, pretype), convert(Cint, maxl))
end

function __SUNSPFGMRSetPrecType(S::SUNLinearSolver, pretype::Cint)
    ccall((:SUNSPFGMRSetPrecType, libsundials_sunlinsolspfgmr), Cint, (SUNLinearSolver, Cint), S, pretype)
end

function SUNSPFGMRSetPrecType(S, pretype)
    __SUNSPFGMRSetPrecType(S, convert(Cint, pretype))
end

function __SUNSPFGMRSetGSType(S::SUNLinearSolver, gstype::Cint)
    ccall((:SUNSPFGMRSetGSType, libsundials_sunlinsolspfgmr), Cint, (SUNLinearSolver, Cint), S, gstype)
end

function SUNSPFGMRSetGSType(S, gstype)
    __SUNSPFGMRSetGSType(S, convert(Cint, gstype))
end

function __SUNSPFGMRSetMaxRestarts(S::SUNLinearSolver, maxrs::Cint)
    ccall((:SUNSPFGMRSetMaxRestarts, libsundials_sunlinsolspfgmr), Cint, (SUNLinearSolver, Cint), S, maxrs)
end

function SUNSPFGMRSetMaxRestarts(S, maxrs)
    __SUNSPFGMRSetMaxRestarts(S, convert(Cint, maxrs))
end

function SUNLinSolGetType_SPFGMR(S::SUNLinearSolver)
    ccall((:SUNLinSolGetType_SPFGMR, libsundials_sunlinsolspfgmr), SUNLinearSolver_Type, (SUNLinearSolver,), S)
end

function SUNLinSolInitialize_SPFGMR(S::SUNLinearSolver)
    ccall((:SUNLinSolInitialize_SPFGMR, libsundials_sunlinsolspfgmr), Cint, (SUNLinearSolver,), S)
end

function SUNLinSolSetATimes_SPFGMR(S::SUNLinearSolver, A_data, ATimes::ATimesFn)
    ccall((:SUNLinSolSetATimes_SPFGMR, libsundials_sunlinsolspfgmr), Cint, (SUNLinearSolver, Ptr{Cvoid}, ATimesFn), S, A_data, ATimes)
end

function SUNLinSolSetPreconditioner_SPFGMR(S::SUNLinearSolver, P_data, Pset::PSetupFn, Psol::PSolveFn)
    ccall((:SUNLinSolSetPreconditioner_SPFGMR, libsundials_sunlinsolspfgmr), Cint, (SUNLinearSolver, Ptr{Cvoid}, PSetupFn, PSolveFn), S, P_data, Pset, Psol)
end

function __SUNLinSolSetScalingVectors_SPFGMR(S::SUNLinearSolver, s1::N_Vector, s2::N_Vector)
    ccall((:SUNLinSolSetScalingVectors_SPFGMR, libsundials_sunlinsolspfgmr), Cint, (SUNLinearSolver, N_Vector, N_Vector), S, s1, s2)
end

function SUNLinSolSetScalingVectors_SPFGMR(S, s1, s2)
    __s1 = convert(NVector, s1)
    __s2 = convert(NVector, s2)
    __SUNLinSolSetScalingVectors_SPFGMR(S, convert(N_Vector, __s1), convert(N_Vector, __s2))
end

function SUNLinSolSetup_SPFGMR(S::SUNLinearSolver, A::SUNMatrix)
    ccall((:SUNLinSolSetup_SPFGMR, libsundials_sunlinsolspfgmr), Cint, (SUNLinearSolver, SUNMatrix), S, A)
end

function __SUNLinSolSolve_SPFGMR(S::SUNLinearSolver, A::SUNMatrix, x::N_Vector, b::N_Vector, tol::realtype)
    ccall((:SUNLinSolSolve_SPFGMR, libsundials_sunlinsolspfgmr), Cint, (SUNLinearSolver, SUNMatrix, N_Vector, N_Vector, realtype), S, A, x, b, tol)
end

function SUNLinSolSolve_SPFGMR(S, A, x, b, tol)
    __x = convert(NVector, x)
    __b = convert(NVector, b)
    __SUNLinSolSolve_SPFGMR(S, A, convert(N_Vector, __x), convert(N_Vector, __b), tol)
end

function SUNLinSolNumIters_SPFGMR(S::SUNLinearSolver)
    ccall((:SUNLinSolNumIters_SPFGMR, libsundials_sunlinsolspfgmr), Cint, (SUNLinearSolver,), S)
end

function SUNLinSolResNorm_SPFGMR(S::SUNLinearSolver)
    ccall((:SUNLinSolResNorm_SPFGMR, libsundials_sunlinsolspfgmr), realtype, (SUNLinearSolver,), S)
end

function SUNLinSolResid_SPFGMR(S::SUNLinearSolver)
    ccall((:SUNLinSolResid_SPFGMR, libsundials_sunlinsolspfgmr), N_Vector, (SUNLinearSolver,), S)
end

function SUNLinSolLastFlag_SPFGMR(S::SUNLinearSolver)
    ccall((:SUNLinSolLastFlag_SPFGMR, libsundials_sunlinsolspfgmr), Clong, (SUNLinearSolver,), S)
end

function SUNLinSolSpace_SPFGMR(S::SUNLinearSolver, lenrwLS, leniwLS)
    ccall((:SUNLinSolSpace_SPFGMR, libsundials_sunlinsolspfgmr), Cint, (SUNLinearSolver, Ptr{Clong}, Ptr{Clong}), S, lenrwLS, leniwLS)
end

function SUNLinSolFree_SPFGMR(S::SUNLinearSolver)
    ccall((:SUNLinSolFree_SPFGMR, libsundials_sunlinsolspfgmr), Cint, (SUNLinearSolver,), S)
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sunlinsol/sunlinsol_spgmr.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __SUNSPGMR(y::N_Vector, pretype::Cint, maxl::Cint)
    ccall((:SUNSPGMR, libsundials_sunlinsolspgmr), SUNLinearSolver, (N_Vector, Cint, Cint), y, pretype, maxl)
end

function SUNSPGMR(y, pretype, maxl)
    __y = convert(NVector, y)
    __SUNSPGMR(convert(N_Vector, __y), convert(Cint, pretype), convert(Cint, maxl))
end

function __SUNSPGMRSetPrecType(S::SUNLinearSolver, pretype::Cint)
    ccall((:SUNSPGMRSetPrecType, libsundials_sunlinsolspgmr), Cint, (SUNLinearSolver, Cint), S, pretype)
end

function SUNSPGMRSetPrecType(S, pretype)
    __SUNSPGMRSetPrecType(S, convert(Cint, pretype))
end

function __SUNSPGMRSetGSType(S::SUNLinearSolver, gstype::Cint)
    ccall((:SUNSPGMRSetGSType, libsundials_sunlinsolspgmr), Cint, (SUNLinearSolver, Cint), S, gstype)
end

function SUNSPGMRSetGSType(S, gstype)
    __SUNSPGMRSetGSType(S, convert(Cint, gstype))
end

function __SUNSPGMRSetMaxRestarts(S::SUNLinearSolver, maxrs::Cint)
    ccall((:SUNSPGMRSetMaxRestarts, libsundials_sunlinsolspgmr), Cint, (SUNLinearSolver, Cint), S, maxrs)
end

function SUNSPGMRSetMaxRestarts(S, maxrs)
    __SUNSPGMRSetMaxRestarts(S, convert(Cint, maxrs))
end

function SUNLinSolGetType_SPGMR(S::SUNLinearSolver)
    ccall((:SUNLinSolGetType_SPGMR, libsundials_sunlinsolspgmr), SUNLinearSolver_Type, (SUNLinearSolver,), S)
end

function SUNLinSolInitialize_SPGMR(S::SUNLinearSolver)
    ccall((:SUNLinSolInitialize_SPGMR, libsundials_sunlinsolspgmr), Cint, (SUNLinearSolver,), S)
end

function SUNLinSolSetATimes_SPGMR(S::SUNLinearSolver, A_data, ATimes::ATimesFn)
    ccall((:SUNLinSolSetATimes_SPGMR, libsundials_sunlinsolspgmr), Cint, (SUNLinearSolver, Ptr{Cvoid}, ATimesFn), S, A_data, ATimes)
end

function SUNLinSolSetPreconditioner_SPGMR(S::SUNLinearSolver, P_data, Pset::PSetupFn, Psol::PSolveFn)
    ccall((:SUNLinSolSetPreconditioner_SPGMR, libsundials_sunlinsolspgmr), Cint, (SUNLinearSolver, Ptr{Cvoid}, PSetupFn, PSolveFn), S, P_data, Pset, Psol)
end

function __SUNLinSolSetScalingVectors_SPGMR(S::SUNLinearSolver, s1::N_Vector, s2::N_Vector)
    ccall((:SUNLinSolSetScalingVectors_SPGMR, libsundials_sunlinsolspgmr), Cint, (SUNLinearSolver, N_Vector, N_Vector), S, s1, s2)
end

function SUNLinSolSetScalingVectors_SPGMR(S, s1, s2)
    __s1 = convert(NVector, s1)
    __s2 = convert(NVector, s2)
    __SUNLinSolSetScalingVectors_SPGMR(S, convert(N_Vector, __s1), convert(N_Vector, __s2))
end

function SUNLinSolSetup_SPGMR(S::SUNLinearSolver, A::SUNMatrix)
    ccall((:SUNLinSolSetup_SPGMR, libsundials_sunlinsolspgmr), Cint, (SUNLinearSolver, SUNMatrix), S, A)
end

function __SUNLinSolSolve_SPGMR(S::SUNLinearSolver, A::SUNMatrix, x::N_Vector, b::N_Vector, tol::realtype)
    ccall((:SUNLinSolSolve_SPGMR, libsundials_sunlinsolspgmr), Cint, (SUNLinearSolver, SUNMatrix, N_Vector, N_Vector, realtype), S, A, x, b, tol)
end

function SUNLinSolSolve_SPGMR(S, A, x, b, tol)
    __x = convert(NVector, x)
    __b = convert(NVector, b)
    __SUNLinSolSolve_SPGMR(S, A, convert(N_Vector, __x), convert(N_Vector, __b), tol)
end

function SUNLinSolNumIters_SPGMR(S::SUNLinearSolver)
    ccall((:SUNLinSolNumIters_SPGMR, libsundials_sunlinsolspgmr), Cint, (SUNLinearSolver,), S)
end

function SUNLinSolResNorm_SPGMR(S::SUNLinearSolver)
    ccall((:SUNLinSolResNorm_SPGMR, libsundials_sunlinsolspgmr), realtype, (SUNLinearSolver,), S)
end

function SUNLinSolResid_SPGMR(S::SUNLinearSolver)
    ccall((:SUNLinSolResid_SPGMR, libsundials_sunlinsolspgmr), N_Vector, (SUNLinearSolver,), S)
end

function SUNLinSolLastFlag_SPGMR(S::SUNLinearSolver)
    ccall((:SUNLinSolLastFlag_SPGMR, libsundials_sunlinsolspgmr), Clong, (SUNLinearSolver,), S)
end

function SUNLinSolSpace_SPGMR(S::SUNLinearSolver, lenrwLS, leniwLS)
    ccall((:SUNLinSolSpace_SPGMR, libsundials_sunlinsolspgmr), Cint, (SUNLinearSolver, Ptr{Clong}, Ptr{Clong}), S, lenrwLS, leniwLS)
end

function SUNLinSolFree_SPGMR(S::SUNLinearSolver)
    ccall((:SUNLinSolFree_SPGMR, libsundials_sunlinsolspgmr), Cint, (SUNLinearSolver,), S)
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sunlinsol/sunlinsol_sptfqmr.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __SUNSPTFQMR(y::N_Vector, pretype::Cint, maxl::Cint)
    ccall((:SUNSPTFQMR, libsundials_sunlinsolsptfqmr), SUNLinearSolver, (N_Vector, Cint, Cint), y, pretype, maxl)
end

function SUNSPTFQMR(y, pretype, maxl)
    __y = convert(NVector, y)
    __SUNSPTFQMR(convert(N_Vector, __y), convert(Cint, pretype), convert(Cint, maxl))
end

function __SUNSPTFQMRSetPrecType(S::SUNLinearSolver, pretype::Cint)
    ccall((:SUNSPTFQMRSetPrecType, libsundials_sunlinsolsptfqmr), Cint, (SUNLinearSolver, Cint), S, pretype)
end

function SUNSPTFQMRSetPrecType(S, pretype)
    __SUNSPTFQMRSetPrecType(S, convert(Cint, pretype))
end

function __SUNSPTFQMRSetMaxl(S::SUNLinearSolver, maxl::Cint)
    ccall((:SUNSPTFQMRSetMaxl, libsundials_sunlinsolsptfqmr), Cint, (SUNLinearSolver, Cint), S, maxl)
end

function SUNSPTFQMRSetMaxl(S, maxl)
    __SUNSPTFQMRSetMaxl(S, convert(Cint, maxl))
end

function SUNLinSolGetType_SPTFQMR(S::SUNLinearSolver)
    ccall((:SUNLinSolGetType_SPTFQMR, libsundials_sunlinsolsptfqmr), SUNLinearSolver_Type, (SUNLinearSolver,), S)
end

function SUNLinSolInitialize_SPTFQMR(S::SUNLinearSolver)
    ccall((:SUNLinSolInitialize_SPTFQMR, libsundials_sunlinsolsptfqmr), Cint, (SUNLinearSolver,), S)
end

function SUNLinSolSetATimes_SPTFQMR(S::SUNLinearSolver, A_data, ATimes::ATimesFn)
    ccall((:SUNLinSolSetATimes_SPTFQMR, libsundials_sunlinsolsptfqmr), Cint, (SUNLinearSolver, Ptr{Cvoid}, ATimesFn), S, A_data, ATimes)
end

function SUNLinSolSetPreconditioner_SPTFQMR(S::SUNLinearSolver, P_data, Pset::PSetupFn, Psol::PSolveFn)
    ccall((:SUNLinSolSetPreconditioner_SPTFQMR, libsundials_sunlinsolsptfqmr), Cint, (SUNLinearSolver, Ptr{Cvoid}, PSetupFn, PSolveFn), S, P_data, Pset, Psol)
end

function __SUNLinSolSetScalingVectors_SPTFQMR(S::SUNLinearSolver, s1::N_Vector, s2::N_Vector)
    ccall((:SUNLinSolSetScalingVectors_SPTFQMR, libsundials_sunlinsolsptfqmr), Cint, (SUNLinearSolver, N_Vector, N_Vector), S, s1, s2)
end

function SUNLinSolSetScalingVectors_SPTFQMR(S, s1, s2)
    __s1 = convert(NVector, s1)
    __s2 = convert(NVector, s2)
    __SUNLinSolSetScalingVectors_SPTFQMR(S, convert(N_Vector, __s1), convert(N_Vector, __s2))
end

function SUNLinSolSetup_SPTFQMR(S::SUNLinearSolver, A::SUNMatrix)
    ccall((:SUNLinSolSetup_SPTFQMR, libsundials_sunlinsolsptfqmr), Cint, (SUNLinearSolver, SUNMatrix), S, A)
end

function __SUNLinSolSolve_SPTFQMR(S::SUNLinearSolver, A::SUNMatrix, x::N_Vector, b::N_Vector, tol::realtype)
    ccall((:SUNLinSolSolve_SPTFQMR, libsundials_sunlinsolsptfqmr), Cint, (SUNLinearSolver, SUNMatrix, N_Vector, N_Vector, realtype), S, A, x, b, tol)
end

function SUNLinSolSolve_SPTFQMR(S, A, x, b, tol)
    __x = convert(NVector, x)
    __b = convert(NVector, b)
    __SUNLinSolSolve_SPTFQMR(S, A, convert(N_Vector, __x), convert(N_Vector, __b), tol)
end

function SUNLinSolNumIters_SPTFQMR(S::SUNLinearSolver)
    ccall((:SUNLinSolNumIters_SPTFQMR, libsundials_sunlinsolsptfqmr), Cint, (SUNLinearSolver,), S)
end

function SUNLinSolResNorm_SPTFQMR(S::SUNLinearSolver)
    ccall((:SUNLinSolResNorm_SPTFQMR, libsundials_sunlinsolsptfqmr), realtype, (SUNLinearSolver,), S)
end

function SUNLinSolResid_SPTFQMR(S::SUNLinearSolver)
    ccall((:SUNLinSolResid_SPTFQMR, libsundials_sunlinsolsptfqmr), N_Vector, (SUNLinearSolver,), S)
end

function SUNLinSolLastFlag_SPTFQMR(S::SUNLinearSolver)
    ccall((:SUNLinSolLastFlag_SPTFQMR, libsundials_sunlinsolsptfqmr), Clong, (SUNLinearSolver,), S)
end

function SUNLinSolSpace_SPTFQMR(S::SUNLinearSolver, lenrwLS, leniwLS)
    ccall((:SUNLinSolSpace_SPTFQMR, libsundials_sunlinsolsptfqmr), Cint, (SUNLinearSolver, Ptr{Clong}, Ptr{Clong}), S, lenrwLS, leniwLS)
end

function SUNLinSolFree_SPTFQMR(S::SUNLinearSolver)
    ccall((:SUNLinSolFree_SPTFQMR, libsundials_sunlinsolsptfqmr), Cint, (SUNLinearSolver,), S)
end


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

using Documenter
using FEMMaterials

deploydocs(
    repo = "github.com/JuliaFEM/FEMMaterials.jl.git",
    target = "build",
    deps = nothing,
    make = nothing)


text:
using JULES.Operators: ∂xᶠᵃᵃ, ∂yᵃᶠᵃ, ∂zᵃᵃᶠ

####
#### Pressure gradient ∇p terms for entropy S = ρs
####

@inline ∂p∂x(i, j, k, grid, tvar, gases, gravity, ρ, ρũ, ρc̃) = ∂xᶠᵃᵃ(i, j, k, grid, diagnose_pressure, tvar, gases, gravity, ρ, ρũ, ρc̃)
@inline ∂p∂y(i, j, k, grid, tvar, gases, gravity, ρ, ρũ, ρc̃) = ∂yᵃᶠᵃ(i, j, k, grid, diagnose_pressure, tvar, gases, gravity, ρ, ρũ, ρc̃)
@inline ∂p∂z(i, j, k, grid, tvar, gases, gravity, ρ, ρũ, ρc̃) = ∂zᵃᵃᶠ(i, j, k, grid, diagnose_pressure, tvar, gases, gravity, ρ, ρũ, ρc̃)


text:
using JuLIP: set_transient!, get_data, has_data

using LinearAlgebra: UniformScaling, eigen, dot, Hermitian, mul!, rmul!

# this file implements the standard spectral decomposition
# calculator for energy, forces, etc.

# ================= diagonalisation ====================
#
# we discussed whether to distinguish k = 0 (Real, Symmetric),
# k ≠ 0 (Complex, Hermitian). For large systems with only Gamma-point,
# this might seem useful, but it actually seems that the timings are
# not so different, e.g., for a 1000 x 1000 system,
#    Float64 ~ 400ms, ComplexF64 ~ 500ms per `eig`.
# this indicates that we probably shouldn't bother.
# This is difference, btw, for LU decompositions, so for the Contour
# calculator we will most likely

eig_(args...) = (eigen(args...)...,)

"""
`sorted_eig`:  helper function to compute eigenvalues, then sort them
in ascending order and sort the eig-vectors as well.
"""
function sorted_eig(H::Matrix, M::Matrix)
   epsn, C = eig_(Hermitian(H), Hermitian(M))
   Isort = sortperm(epsn)
   return epsn[Isort], C[:, Isort]
end

function sorted_eig(H::Matrix, ::UniformScaling)
   epsn, C = eig_(Hermitian(H))
   Isort = sortperm(epsn)
   return epsn[Isort], C[:, Isort]
end

"""
`update_eig!(atm::AbstractAtoms, tbm::TBModel)` : updates the hamiltonians
and spectral decompositions on the MP grid.
"""
function update_eig!(atm::AbstractAtoms, H::SparseSKH{ISORTH}, tbm::TBModel) where {ISORTH}
   wrk = _alloc_full(H)
   for (w, k) in tbm.bzquad
      if !has_k_array(atm, :epsn, k)
         Hf, Mf = full!(wrk, H, k)
         epsn, C = sorted_eig(Hf, Mf)
         # TODO: probably remove this - but need to rewrite site_energy first
         set_k_array!(atm, Mf, :M, k)
         set_k_array!(atm, epsn, :epsn, k)
         set_k_array!(atm, C, :C, k)
      end
   end
end

# ==================== update all arrays =================
#    e.g., prior to an energy or force calculation

"""
`update!(atm::AbstractAtoms, tbm:TBModel)`: checks whether the precomputed
data stored in `at` is still there - JuLIP deletes all arrays when
atom positions are updated.

* spectral decompositions (`update_eig!`)
* the fermi-level (`update_eF!`)
"""
function update!(at::AbstractAtoms, tbm::TBModel)
   # check whether the :tbupdateflag exists;  if it does then atom positions
   # have not changed since last time that update! was called
   if has_data(at, :tbupdateflag)
      return nothing
   end
   # set the update flag (will be deleted as soon as atom positions change)
   # ( we do this *before* the update so that we don't go into an infinite
   #    loop with the update trying to update; see e.g. BZiter )
   set_transient!(at, :tbupdateflag, 0)
   # if the flag does not exist, then we update everything
   skh = SparseSKH(tbm.H, at)  # this also stores skh for later use
   update_eig!(at, skh, tbm)
   update!(at, tbm.potential, tbm)
   return nothing
end

# ================ Density Matrix and Energy ================

"""
`densitymatrix(tbm, at) -> Γ`:

### Input

* `tbm::TBModel` : calculator
* `at::AbstractAtoms` : configuration

### Output

`Γ::Matrix{Float64}`: density matrix Γ = Σ_k w^k Σ_s f(ϵ_s^k) ψ_s^k ⊗ ψ_s^k
and f(ϵ) is the occupancy
"""
function densitymatrix(tbm::TBModel, at::AbstractAtoms)
   N = ndofs(tbm.H, at)
   Γ = zeros(N, N)
   for (w, _, ϵ, ψ) in BZiter(tbm, at)
      fs = occupancy(tbm.potential, ϵ)
      for a = 1:N, b = 1:N
         Γ[a,b] += w * fs * real(ψ[a] * ψ[b]')
      end
   end
   return Γ
end

# this is imported from JuLIP
_energy(tbm::TBModel, at::AbstractAtoms) = (
   sum( w * energy(tbm.potential, ϵ) for (w, _1, ϵ, _2) in BZiter(tbm, at) )
   + energy(tbm.Vrep, at) )

# function energy(tbm::TBModel, at::AbstractAtoms)
#    update!(at, tbm)
#    Etot = 0.0
#    for (w, k) in tbm.bzquad
#       epsn_k = get_k_array(at, :epsn, k)
#       Etot += w * sum( energy(tbm.potential, epsn_k) )
#    end
#    return Etot # + energy(tbm.Vrep, at)
# end

# ========================== Forces ==========================

# this is an old force computation that requires a SKHamiltonian structure
#    F_n = - ∑_s f'(ϵ_s) < ψ_s | H,n - ϵ_s * M,n | ψ_s >
# but instead of computing H,n, M,n as matrices, this assembly loops over
# the non-zero blocks first and the inner loop is over the derivative ,n
# it prevents allocating lots of sparse matrices for the hamiltonian
# derivatives.

function _forces_k!(frc::AbstractVector{JVecF},
                                 at::AbstractAtoms, tbm::TBModel,
                                 H::SKHamiltonian{ISORTH,NORB}, k::JVecF,
                                 skhg, w) where {ISORTH, NORB}
   # obtain the precomputed arrays
   epsn = get_k_array(at, :epsn, k)::Vector{Float64}
   C = get_k_array(at, :C, k)::Matrix{ComplexF64}
   df = grad(tbm.potential, epsn)::Vector{Float64}

   # precompute some products
   # TODO: optimise these two lines? These are O(N^3) scaling, while overall
   #       force assembly should be just O(N^2) (but can we beat BLAS3?)
   C_df_Ct = (C * (df' .* C)')
   C_dfepsn_Ct = (C * ((df.*epsn)' .* C)')

   # an array replacing dM_ij when the model is orthogonal
   dM_ij = zero(typeof(skhg.dH[1]))

   for n = 1:length(skhg.i)
      i, j, dH_ij, dH_ii, S = skhg.i[n], skhg.j[n], skhg.dH[n], skhg.dOS[n], skhg.Rcell[n]
      if !ISORTH; dM_ij = skhg.dM[n]; end
      Ii, Ij = indexblock(i, H), indexblock(j, H)
      eikr = exp(im * dot(S, k))::Complex{Float64}

      @inbounds for a = 1:NORB, b = 1:NORB
         t1 = 2.0 * real(C_df_Ct[Ii[a], Ij[b]] * eikr)
         t2 = 2.0 * real(C_dfepsn_Ct[Ii[a], Ij[b]] * eikr)
         t3 = real(C_df_Ct[Ii[a],Ii[b]])
         frc[j] += (-t1) * dH_ij[:,a,b] + t2 * dM_ij[:,a,b] - t3 * dH_ii[:,a,b]
         frc[i] += t3 * dH_ii[:,a,b]
      end
   end

   return frc
end

# * `forces` is imported from JuLIP
# * this implementation is the old version from Atoms.jl, which makes
#   specific assumptions about the structure of the hamiltonian, hence is
#   only valid for SK-type hamiltonians.
#
function _forces(tbm::TBModel{HT}, atm::AbstractAtoms) where {HT <: SKHamiltonian}
   update!(atm, tbm)
   skhg = SparseSKHgrad(tbm.H, atm)
   frc = forces(tbm.Vrep, atm)
   for (w, k) in tbm.bzquad
      _forces_k!(frc, atm, tbm, tbm.H, k, skhg, w)
   end
   return frc
end

function energy(tbm::TBModel,
                at::AbstractAtoms;
                domain = 1:length(at))
   if domain == 1:length(at); return _energy(tbm, at); end
   update!(at, tbm)
   In0 = indexblock(domain, tbm.H)
   E = 0.0
   for (w, k) in tbm.bzquad
      epsn_k = get_k_array(at, :epsn, k)
      M_k = get_k_array(at, :M, k)
      C_k = get_k_array(at, :C, k)
      MC_k = isorth(tbm) ? C_k[In0,:] : (M_k[In0, :] * C_k)
      ψ² = sum( conj(C_k[In0, :] .* MC_k), dims=1 )[:]
      E += w * sum(energy(tbm.potential, epsn_k)  .* ψ²)
   end
   return real(E) + energy(tbm.Vrep, at; domain=domain)
end

site_energy(tbm::TBModel, at::AbstractAtoms, n0::Integer) =
      energy(tbm, at; domain = [n0])

# ================ Band-structure =================

"""
`band_structure(tbm::TBModel, at::AbstractAtoms) -> k, E`

where
* `k` is a list of k-points
* `E` is a ndofs x nk matrix of energies
"""
function band_structure(tbm::TBModel, at::AbstractAtoms)
   update!(at, tbm)
   na = ndofs(tbm.H, at)
   K = JVecF[]
   E = Float64[]
   for (w, k) in tbm.bzquad
      push!(K, k)
      append!(E, get_k_array(at, :epsn, k))
   end
   return K, reshape(E, na, length(E) ÷ na)
end

"""
`spectrum(tbm, at) -> E::Vector`
"""
spectrum(tbm, at) = band_structure(tbm, at)[2][:]

function spectrum(k::JVec, tbm::TBModel, at::AbstractAtoms)
   H = SparseSKH(tbm.H, at)
   wrk = _alloc_full(H)
   Hk, Mk = full!(wrk, H, k)
   epsn, C = sorted_eig(Hk, Mk)
   return epsn
end

function band_structure(K::AbstractVector, tbm::TBModel, at::AbstractAtoms)
   EE = [ spectrum(k, tbm, at) for k in K ]
   bands = [ zeros(length(EE)) for n = 1:length(EE[1]) ]
   for n = 1:length(EE[1]), m = 1:length(EE)
      bands[n][m] = EE[m][n]
   end
   return bands
end

# """
# `band_structure_near_eF(Nb::Integer, at::AbstractAtoms, tbm::TBModel) -> k, E`
#
# get 2*Nb+1 bands around the fermi level
# """
# function band_structure_near_eF(Nb, at::AbstractAtoms, tbm::TBModel)
#    update!(at, tbm)
#    Ne = ndofs(tbm.H, at)
#    K = JVecF[]
#    E = Float64[]
#    eF = get_eF(tbm.potential)
#    for (w, k) in tbm.bzquad
#       push!(K, k)
#       epsn_k = get_k_array(tbm, :epsn, k)
#
#       E[Nb+1,n] = epsn_k[nf]
#       for j = 1:Nb
#          E[Nb+1-j,n] = epsn_k[nf-j]
#          E[Nb+1+j,n] = epsn_k[nf+j]
#       end
#    end
#    return K, E
# end

# =================== Site Forces for a given k-point =======================
# The derivatives of the site energy is computed by dual technique.

function _dEs_k!(dEs::AbstractVector{JVecF},
                  at::AbstractAtoms, tbm::TBModel,
                  H::SKHamiltonian{ISORTH,NORB}, k::JVecF,
                  skhg, w, Idom::AbstractVector{TI}
         ) where {ISORTH, NORB, TI <: Integer}

   # obtain the orbital-indices corresponding to the site indices `Idom`
   In0  = indexblock(Idom, tbm.H)
   # obtain the precomputed arrays
   epsn = get_k_array(at, :epsn, k)::Vector{Float64}
   C    = get_k_array(at, :C, k)::Matrix{ComplexF64}
   f    = energy(tbm.potential, epsn)::Vector{Float64}
   df   = grad(tbm.potential, epsn)::Vector{Float64}
   M    = get_k_array(at, :M, k)::Matrix{ComplexF64}
   MC   = isorth(tbm) ? C[In0,:] : (M[In0, :] * C)
   ψ²   = sum( conj(C[In0, :]) .* MC , dims=1 )[:]
   # precompute some products
   C_f_Ct         = (C * (f' .* C)')
   C_f_ψ²_Ct      = (C * ( (f .* ψ²)' .* C )')
   C_df_ψ²_Ct     = (C * ( (df .* ψ²)' .* C )')
   C_dfepsn_ψ²_Ct = (C * ( (df .* epsn .* ψ²)' .* C )')
   # an array replacing dM_ij when the model is orthogonal
   dM_ij = zero(typeof(skhg.dH[1]))
   _A1 = zeros(ComplexF64, size(C))

   # precompute pseudo-inverse:
   Nelc = length(epsn)
   diff_f_inv_ψst = zeros(ComplexF64, Nelc, Nelc)
   diff_fepsn_inv_ψst = zeros(ComplexF64, Nelc, Nelc)
	for t = 1:Nelc, s = 1:Nelc
      # _aa = sum( conj(C[In0[ii], s]) * MC[ii, t] for ii = 1:length(In0))
      _aa = zero(ComplexF64)
      for ii = 1:length(In0)
         _aa += conj(C[In0[ii], s]) * MC[ii, t]
      end

		if abs(epsn[t]-epsn[s]) > 1e-10
        	diff_f_inv_ψst[t,s] =
               (f[s]-f[t]) * _aa/ (epsn[s]-epsn[t])
         diff_fepsn_inv_ψst[t,s] =
               (f[s]*epsn[s]-f[t]*epsn[t]) * _aa / (epsn[s]-epsn[t])
         # diff_f_inv_ψst[t,s] =
         #       (f[s]-f[t]) * ( C[In0,s]' * MC[:,t] )[1] / (epsn[s]-epsn[t])
         # diff_fepsn_inv_ψst[t,s] =
         #       (f[s]*epsn[s]-f[t]*epsn[t]) * ( C[In0,s]' * MC[:,t] )[1] / (epsn[s]-epsn[t])
      else
        	diff_f_inv_ψst[t,s] = df[s] * _aa
         diff_fepsn_inv_ψst[t,s] = (df[s]*epsn[s] + f[s]) * _aa
         # diff_f_inv_ψst[t,s] = df[s] * ( C[In0,s]' * MC[:,t] )[1]
         # diff_fepsn_inv_ψst[t,s] = (df[s]*epsn[s] + f[s]) * ( C[In0,s]' * MC[:,t] )[1]
      end
	end
   # precompute the arrays (H-ϵ_s⋅M)⁻[ψ_s]
   # pinvC_f_Ct = (C * diff_f_inv_ψst * C')
   pinvC_f_Ct = mul!(_A1, C, diff_f_inv_ψst) * adjoint(C)
   # pinvC_fepsn_Ct = (C * diff_fepsn_inv_ψst * C')
   pinvC_fepsn_Ct = mul!(_A1, C, diff_fepsn_inv_ψst) * adjoint(C)

   for n = 1:length(skhg.i)
      i, j, dH_ij, dH_ii, S = skhg.i[n], skhg.j[n], skhg.dH[n], skhg.dOS[n], skhg.Rcell[n]
      if !ISORTH; dM_ij = skhg.dM[n]; end
      Ii, Ij = indexblock(i, H), indexblock(j, H)
      eikr = exp(im * dot(S, k))::Complex{Float64}
      @inbounds for a = 1:NORB, b = 1:NORB
         if i ∈ Idom
            t4  = real(C_f_Ct[Ii[a], Ij[b]] * eikr)
            dEs[j] += t4 * dM_ij[:,a,b]
            dEs[i] += - t4 * dM_ij[:,a,b]
         end
         # be careful: the matrix pinvC_f_Ct is not symmetric
         t5 = real(pinvC_f_Ct[Ii[a],Ij[b]] * eikr) +
              real(pinvC_f_Ct[Ij[b],Ii[a]] * eikr)
         t6 = real(pinvC_fepsn_Ct[Ii[a],Ij[b]] * eikr) +
              real(pinvC_fepsn_Ct[Ij[b],Ii[a]] * eikr)
         t7 = real(pinvC_f_Ct[Ii[a],Ii[b]])
         dEs[j] += t5 * dH_ij[:,a,b] - t6 * dM_ij[:,a,b] + t7 * dH_ii[:,a,b]
         dEs[i] += - t7 * dH_ii[:,a,b]
      end
   end
   return dEs
end

# Derivative of site energy using dual technique
#
function forces(tbm::TBModel{HT}, atm::AbstractAtoms{T};
                domain = 1:length(atm) ) where {T, HT <: SKHamiltonian}
   if domain ==  1:length(atm); return _forces(tbm, atm); end
   update!(atm, tbm)
   skhg = SparseSKHgrad(tbm.H, atm)
   dE = zeros(JVec{T}, length(atm))
   for (w, k) in tbm.bzquad
      _dEs_k!(dE, atm, tbm, tbm.H, k, skhg, w, domain)
   end
   return rmul!(dE + forces(tbm.Vrep, atm; domain  = domain), -1)
end

site_energy_d(tbm::TBModel{HT}, atm::AbstractAtoms,
               n0::Integer) where {HT <: SKHamiltonian} =
   rmul!(forces(tbm, atm; domain = [n0]), -1)

# ========================== Virial ==========================

# compute the virial ( = - ∂E / ∂F ) through
# vir = - ∑_s f'(ϵ_s) < ψ_s | ∂H/∂u_{αβ} - ϵ_s * ∂M/∂u_{αβ} | ψ_s >
# but instead of computing H and M as matrices, this assembly loops over
# the non-zero blocks first and the inner loop is over the derivative

# TODO: add domain dependence to virial calculations. (maybe)

function _virial_k!(vir_k::Array{Float64,2},
                           at::AbstractAtoms, tbm::TBModel,
                           H::SKHamiltonian{ISORTH,NORB}, k::JVecF,
                           skhg, w) where {ISORTH, NORB}
   # deformation
   F = cell(at)'
   F⁻ = inv(F)
   X = positions(at)|>mat

   # obtain the precomputed arrays
   epsn = get_k_array(at, :epsn, k)::Vector{Float64}
   C = get_k_array(at, :C, k)::Matrix{ComplexF64}
   df = grad(tbm.potential, epsn)::Vector{Float64}

   # precompute some products
   # TODO: optimise these two lines? These are O(N^3) scaling, while overall
   #       force assembly should be just O(N^2) (but can we beat BLAS3?)
   C_df_Ct = (C * (df' .* C)')
   C_dfepsn_Ct = (C * ((df.*epsn)' .* C)')

   # an array replacing dM_ij when the model is orthogonal
   dM_ij = zero(typeof(skhg.dH[1]))
   dH_ij_F = zeros(Float64, 3, size(skhg.dH[1])[1], size(skhg.dH[1])[2], size(skhg.dH[1])[3])
   dH_ii_F = zeros(Float64, 3, size(skhg.dH[1])[1], size(skhg.dH[1])[2], size(skhg.dH[1])[3])
   dM_ij_F = zeros(Float64, 3, size(skhg.dH[1])[1], size(skhg.dH[1])[2], size(skhg.dH[1])[3])

   for n = 1:length(skhg.i)
      i, j, dH_ij, dH_ii, S = skhg.i[n], skhg.j[n], skhg.dH[n], skhg.dOS[n], skhg.Rcell[n]
      if !ISORTH; dM_ij = skhg.dM[n]; end
      Ii, Ij = indexblock(i, H), indexblock(j, H)
      eikr = exp(im * dot(S, k))::Complex{Float64}
      # x ↦ F⋅F⁻x
      F⁻xi = F⁻ * X[:,i]
      F⁻xj = F⁻ * (S + X[:,j])
      for t = 1:3
         dM_ij_F[t,:,:,:] = dM_ij .* ( F⁻xj[t] - F⁻xi[t] )
         dH_ij_F[t,:,:,:] = dH_ij .* ( F⁻xj[t] - F⁻xi[t] )
         dH_ii_F[t,:,:,:] = dH_ii .* ( F⁻xj[t] - F⁻xi[t] )
      end

      @inbounds for a = 1:NORB, b = 1:NORB
         t1 = real(C_df_Ct[Ii[a], Ij[b]] * eikr)
         t2 = real(C_dfepsn_Ct[Ii[a], Ij[b]] * eikr)
         t3 = real(C_df_Ct[Ii[a],Ii[b]])
         vir_k += (-t1) * dH_ij_F[:,:,a,b] + t2 * dM_ij_F[:,:,a,b] - t3 * dH_ii_F[:,:,a,b]
      end
   end

   return vir_k
end

# * `virial` is imported from JuLIP  "CHECK"
# hence is only valid for SK-type hamiltonians.

function virial(tbm::TBModel{HT}, atm::AbstractAtoms) where {HT <: SKHamiltonian}
   update!(atm, tbm)
   skhg = SparseSKHgrad(tbm.H, atm)
   # TODO: fix -- vir = virial(tbm.Vrep, atm)
   vir = zeros(Float64, 3, 3)
   for (w, k) in tbm.bzquad
      vir_k = zeros(Float64, 3, 3)
      vir += w * _virial_k!(vir_k, atm, tbm, tbm.H, k, skhg, w)
   end
   return vir
end


text:
import DataFrames, Statistics

@test standardize([1, 2, 3]) == [-1.0, 0, 1]
@test standardize([1, 2, 3], robust=true) ≈ [-0.67, 0, 0.67] atol=0.01
@test standardize!([1, 2, 3]) == [-1.0, 0, 1]
@test standardize(DataFrames.DataFrame([[1, 2, 3]])) == DataFrames.DataFrame([[-1.0, 0, 1]])
@test size(standardize(DataFrames.DataFrame([[1, missing, 2, 3]]))) == (4, 1)

@test Statistics.mean(perfectNormal(50)) ≈ 0 atol=0.01
@test Statistics.mean(perfectNormal(n = 50)) ≈ 0 atol=0.01


text:
using Documenter, TORA

makedocs(
    modules = [TORA],
    format = Documenter.HTML(),
    sitename = "TORA.jl",
    authors = "Henrique Ferrolho",
    pages = [
        "Home" => "index.md",
        "Getting Started" => [
            "install.md",
            "tutorial.md",
        ],
        "Manual" => [
            "robot.md",
            "problem.md",
            "solver_interfaces.md",
            "solver_log.md",
        ],
        "Advanced" => [
            "new_robot.md",
            "how_it_works.md",
        ]
    ]
)

deploydocs(
    repo = "github.com/JuliaRobotics/TORA.jl",
    devbranch = "main",
)


text:
module species

# TODO: Create custom constructor that creates v_t as part of type

export SpeciesData

"""
A container for species data. The following example creates a species of electrons in SI units.

 * `SpeicesData(mass=me,charge=-el,density=1.0e20,temp=1e3*el)`

SI units are not necessary, but is expected when passing to other methods in PlasmaTools.

"""
type SpeciesData
   mass::Float64
   charge::Float64
   density::Float64
   temperature::Float64
end

# Define constructor so that all keywords must be specified instead of ordered
SpeciesData(;mass=-1.0,charge=-1.0,density=-1.0,temperature=-1.0) =SpeciesData(mass,charge,density,temperature)

end

max_stars_repo_name: gjwilkie/PlasmaTools

max_issues_repo_name: gjwilkie/PlasmaTools

max_forks_repo_name: gjwilkie/PlasmaTools


text:
# Output functions for GeMM

"""
    printheader(io, sep)

Print a list of property names to the given IO stream. This is a helper function
for `dumpinds`.
"""
function printheader(io::IO = stdout, sep::String = "\t")
    #XXX Transfer to a dynamic system? (As in createworld()?)
    print(io, "time", sep)
    print(io, "patch_no", sep)
    print(io, "xloc", sep)
    print(io, "yloc", sep)
    print(io, "temp", sep)
    print(io, "capacity", sep)
    setting("nniches") > 1 && print(io, "prec", sep)
    setting("nniches") > 2 && print(io, "nicheb", sep)
    print(io, "island", sep)
    print(io, "isolation", sep)
    print(io, "invasible", sep)
    print(io, "initpop", sep)
    print(io, "id", sep)
    print(io, "lineage", sep)
    print(io, "new", sep)
    print(io, "tempadaptation", sep)
    print(io, "precadaptation", sep)
    print(io, "size", sep)
    print(io, "lnkgunits", sep)
    print(io, "ngenes", sep)
    for key in setting("traitnames")
        print(io, key, sep)
    end
    println(io)
end

"""
    dumpinds(world, io, sep)

Output all data of individuals in `world` as table to `io`. Columns are separated by `sep`.
WARNING: this produces *very* large files!
"""
function dumpinds(world::Array{Patch, 1}, timestep::Int, io::IO = stdout, sep::String = "\t")
    timestep == 0 && printheader(io, sep)
    onlyisland = setting("static") && timestep > 1
    for patch in world
        #XXX Sometimes, this only dumps three or four individuals, with a population of >10⁴!
        # (Should be fixed)
        (onlyisland && !patch.isisland) && continue
        if !patch.isisland && setting("static")
            community = patch.seedbank
        else
            community = patch.community
        end
        for ind in community
            print(io, timestep, sep)
            print(io, patch.id, sep)
            print(io, patch.location[1], sep)
            print(io, patch.location[2], sep)
            print(io, patch.temp, sep)
            print(io, patch.capacity, sep)
            setting("nniches") > 1 && print(io, patch.prec, sep)
            setting("nniches") > 2 && print(io, patch.nicheb, sep)
            patch.isisland ? print(io, 1, sep) : print(io, 0, sep)
            patch.isolated ? print(io, 1, sep) : print(io, 0, sep)
            patch.invasible ? print(io, 1, sep) : print(io, 0, sep)
            patch.initpop ? print(io, 1, sep) : print(io, 0, sep)
            print(io, ind.id, sep)
            print(io, ind.lineage, sep)
            ind.marked ? print(io, 1, sep) : print(io, 0, sep)
            print(io, ind.tempadaptation, sep)
            print(io, ind.precadaptation, sep)
            print(io, ind.size, sep)
            print(io, length(ind.genome), sep)
            print(io, sum(map(x -> length(x.genes), ind.genome)), sep)
            for key in setting("traitnames")
                try
                    print(io, ind.traits[key], sep)
                catch
                    print(io, "NA", sep)
                end
            end
            println(io)
        end
    end
end

"""
     printheaderid(io, sep)
helper function for indcoord(world, io, sep) to print the column headers.
"""
function printheaderid(io::IO = stdout, sep::String = "\t")
    print(io, "time", sep)
    print(io, "patch_no", sep)
    print(io, "xloc", sep)
    print(io, "yloc", sep)
    print(io, "id", sep)
    println(io)
end

"""
    indcoord(world, io, sep)

Output id and coordinates of individuals in `world` as table to `io`. Columns are separated by `sep`.
"""
function indcoord(world::Array{Patch, 1}, timestep::Int, io::IO = stdout, sep::String = "\t")
    timestep == 0 && printheaderid(io, sep)
    for patch in world
        for ind in patch.community
            print(io, timestep, sep)
            print(io, patch.id, sep)
            print(io, patch.location[1], sep)
            print(io, patch.location[2], sep)
            print(io, ind.id, sep)
            println(io)
        end
    end
end

"""
    makefasta(world, io, onlyisland, sep)

Record the genome of every individual currently alive to the given IO stream.
(High-detail data recording function.) WARNING: this produces *very* large files!
"""
function makefasta(world::Array{Patch, 1}, io::IO = stdout, onlyisland::Bool = false, sep::String = "_")
    for patch in world
        (onlyisland && !patch.isisland) && continue
        if !patch.isisland && setting("static")
            community = patch.seedbank
        else
            community = patch.community
        end
        lineage = ""
        for ind in community
            (!patch.isisland && setting("static") && ind.lineage == lineage) && continue # only one individual per species on mainland
            chrmno = 0
            for chrm in ind.genome
                chrmno += 1
                geneno = 0
                for gene in chrm.genes
                    geneno += 1
                    traits = ""
                    if length(gene.codes) == 0
                        traits *= ","
                    else
                        for trait in gene.codes
                            traits *= string(setting("traitnames")[trait.nameindex]) * ":" * string(trait.value) * ","
                        end
                    end
                    header = ">" * ind.lineage * sep * string(ind.id) * sep * string(chrmno) * sep * string(geneno) * sep * traits
                    if (setting("fasta") == "compat" && occursin("compat", header)) || setting("fasta") != "compat"
                        println(io, header)
                        if setting("compressgenes")
                            println(io, num2seq(gene.sequence))
                        else
                            println(io, gene.sequence)
                        end
                    end
                end
            end
            lineage = ind.lineage
        end
    end
end

"""
    setupdatadir(dir)

Creates the output directory and copies relevant files into it.
"""
function setupdatadir()
    if isdir(setting("dest"))
        #XXX a more useful solution might be to rename "dest" to avoid conflict
        simlog(string(setting("dest"))*" exists. Aborting to avoid overwriting files.", 'e')
    else
        mkpath(setting("dest"))
    end
    simlog("Setting up output directory "*string(setting("dest")))
    writesettings()
    if in("maps", settingkeys())
        for m in setting("maps")
            isempty(m) && continue
            !(isfile(m)) && simlog("Map file "*string(m)*" doesn't exist!", 'e')
            cp(m, joinpath(setting("dest"), basename(m)), force = true)
        end
    end
end

"""
    writesettings()

Record the settings actually used for a simulation run (cf. `getsettings`).
Creates a config file that can be used for future replicate runs.
Also records a time stamp and the current git commit.
"""
function writesettings()
    if isempty(basename(setting("config")))
        settingspath = "settings.conf"
    else
        settingspath = basename(setting("config"))
    end
    open(joinpath(setting("dest"), settingspath), "w") do f
        println(f, "#\n# --- GeMM configuration parameters ---")
        println(f, "# This file was generated automatically.")
        println(f, "# Simulation run on "*string(Dates.format(Dates.now(), "d u Y HH:MM:SS")))
        println(f, "# "*string(split(read(pipeline(`git log`, `head -1`), String), "\n")[1])*"\n")
        for k in settingkeys()
            value = setting(k)
            if isa(value, String)
                value = "\"" * value * "\""
            elseif isa(value, Array)
                vstr = "\""
                        for x in value
                            vstr *= string(x) * ","
                        end
                        value = vstr[1:end-1] * "\""
            end
            println(f, string(k), " ", string(value))
        end
    end
end

"""
    writedata(world, timestep)

Writes simulation output from `world` to table files. (Which
data is recorded depends on the settings.) `timestep` and `setting` information
is used for file name creation.
"""
function writedata(world::Array{Patch,1}, timestep::Int)
    if setting("lineages")
        recordlineages(world, timestep)
        recordstatistics(world)
    end
    if setting("raw")
        filename = "inds_s" * string(setting("seed"))
        filename = joinpath(setting("dest"), filename)
        filename = filename * ".tsv"
        simlog("Writing data \""*string(filename)*"\"")
        open(filename, "a") do file
            dumpinds(world, timestep, file)
        end
    end
    if setting("stats")
        # why not `filename = joinpath(setting("dest"), "pops.tsv")`?
        filename = "pops_s" * string(setting("seed"))
        filename = joinpath(setting("dest"), "pops.tsv")
        simlog("Writing stats to \""*string(filename)*"\"")
        open(filename, "a") do file
            printpopstats(file, world, timestep)
        end
    end
end

"""
    writefasta(world, timestep)

Writes simulation output from `world` to fasta files and an additional file with
just the id and the coordinates of each individual. (Which
data is recorded depends on the settings.) `timestep` and `setting` information
is used for file name creation.
"""
function writefasta(world::Array{Patch,1}, timestep::Int)
    filename = "seqs_s" * string(setting("seed"))
    filename = joinpath(setting("dest"), filename)
    filename = filename * ".fa"
    simlog("Writing fasta \""*string(filename)*"\"")
    open(filename, "a") do file
        makefasta(world, file)
    end
    if setting("dumpindforfasta")
        filename_ind = "inds_s"* string(setting("seed"))
        filename_ind = joinpath(setting("dest"), filename_ind)
        filename_ind = filename_ind * ".tsv"
        simlog("Writing individual data \""*string(filename_ind)*"\"")
        open(filename_ind, "a") do file
            dumpinds(world, timestep, file)
        end
    end
end

"""
    recordstatistics(w)

Write out world properties to the log file for later analysis.
"""
function recordstatistics(world::Array{Patch,1})
    if !isfile(joinpath(setting("dest"), "diversity.log"))
        simlog("population,freespace,lineages,alpha,beta,gamma", 'i', "diversity.log", true)
    end
    popsize = sum(x -> length(x.community), world)
    lineages = unique(reduce(vcat, map(p -> map(x -> x.lineage, p.community), world)))
    div = round.(diversity(world), digits = 3)
    space = freespace(world)
    simlog("Metacommunity size: "*string(popsize)*", lineages: "*string(length(lineages))) # to stdout
    simlog(string(popsize)*","*string(space)*","*string(length(lineages))*","*
           string(div[1])*","*string(div[2])*","*string(div[3]),
           'i', "diversity.log", true)
end

"""
    recordlineages(w)

Save the abundance of each lineage per patch. (Low-detail data recording function.)
"""
function recordlineages(world::Array{Patch,1}, timestep::Int)
    #XXX despite being low-detail, calling this frequently still means a lot of I/O
    if !isfile(joinpath(setting("dest"), "lineages.log"))
        simlog("t,X,Y,lineage,abundance,temp,prec", 'i', "lineages.log", true)
    end
    datastring = ""
    for p in world
        for l in unique(map(x -> x.lineage, p.community))
            datastring *= string(timestep)*","*string(p.location[1])*","*string(p.location[2])*
            ","*string(l)*","*string(length(findall(x -> x.lineage == l, p.community)))*","*
            string(p.temp)*","*string(p.prec)*"\n"
            #TODO the trailing linefeed causes empty lines in the output file between updates
        end
    end
    simlog(datastring, 'i', "lineages.log", true)
end

"""
    printpopheader(io)

Print a list of property names to the given IO stream. This is a helper function
for `printpopstats`.
"""
function printpopheader(io::IO)
    print(io, "time", "\tx", "\ty", "\ttemp", "\tprec", "\tcapacity", "\tisisland",
          "\tlineage", "\tjuveniles", "\tadults", "\ttempadaptationmean",
          "\tprecadaptationmean", "\theterozygosity")
    #XXX this probably should't be hardcoded, use `setting("traitnames")` instead
    traitnames =  ["compat", "compatsd", "dispmean", "dispmeansd", "dispshape", "dispshapesd",
                   "ngenes", "nlnkgunits", "precopt", "precoptsd", "prectol", "prectolsd",
                   "repsize", "repsizesd", "seqsimilarity", "seqsimilaritysd", "seedsize", "seedsizesd",
                   "tempopt", "tempoptsd", "temptol", "temptolsd"]
    for traitname in traitnames
        print(io, "\t", traitname, "mean")
        print(io, "\t", traitname, "std")
    end
    print(io, "\treplicate", "\tconf")
    println(io)
end

"""
    printpopstats(io, world, timestep)

Record statistical information (maximum, minimum, median, standard deviation)
for a range of individual properties, as seen over the whole world population.
(Medium-detail data recording function.)
"""
function printpopstats(io::IO, world::Array{Patch, 1}, timestep::Integer)
    timestep == 0 && printpopheader(io)
    #XXX this probably should't be hardcoded, use `setting("traitnames")` instead
    traitnames =  ["compat", "compatsd", "dispmean", "dispmeansd", "dispshape", "dispshapesd",
                   "ngenes", "nlnkgunits", "precopt", "precoptsd", "prectol", "prectolsd",
                   "repsize", "repsizesd", "seqsimilarity", "seqsimilaritysd", "seedsize", "seedsizesd",
                   "tempopt", "tempoptsd", "temptol", "temptolsd"]
    for patch in world
        lineages = unique(map(i -> i.lineage, patch.community))
        for lineage in lineages
            print(io, timestep, "\t", patch.location[1], "\t", patch.location[2],
                  "\t", patch.temp, "\t", patch.prec, "\t", patch.capacity, "\t", patch.isisland)
            popidxs = findall(i -> i.lineage == lineage, patch.community)
            population = patch.community[popidxs]
            adultidxs = findall(i -> i.size >= i.traits["repsize"], patch.community[popidxs])
            print(io, "\t", population[1].lineage, "\t", length(popidxs) - length(adultidxs), "\t", length(adultidxs),
                  "\t", mean(skipmissing(map(i -> i.tempadaptation, population))),
                  "\t", mean(skipmissing(map(i -> i.precadaptation, population))),
                  "\t", heterozygosity(population))
            for traitname in traitnames
                poptrait = map(i -> i.traits[traitname], population)
                print(io, "\t", mean(skipmissing(poptrait)))
                print(io, "\t", std(skipmissing(poptrait))) # CAVEAT: this returns NaN if only one individual
            end
            print(io, "\t", setting("seed"), "\t", setting("config"))
            println(io)
        end
    end
end

"""
    simlog(msg, category, logfile, onlylog)

Write a log message to STDOUT/STDERR and the specified logfile
(if logging is turned on in the settings).

Categories: `d` (debug), `i` (information, default), `w` (warn), `e` (error)

If `logfile` is the empty string (default: "simulation.log"), the message will
only be printed to the screen. If `onlylog` is true (default: false), the
message is not printed to screen but only to the log.
"""
function simlog(msg::String, category::Char='i', logfile::String="simulation.log", onlylog::Bool=false)
    #TODO Julia now has inbuilt logging facilities: https://docs.julialang.org/en/v1/stdlib/Logging/
    # This function ought to be rewritten to make use of these (especially warning and error macros)
    (isa(category, String) && length(category) == 1) && (category = category[1])
    function logprint(msg::String, tostderr=false)
        if tostderr || !(setting("quiet") || onlylog)
            tostderr ? iostr = stderr : iostr = stdout
            #TODO use `printstyled` to give coloured errors/warnings
            println(iostr, msg)
        end
        if setting("logging")
            #XXX always opening new connections is expensive -> should we buffer the output?
            open(joinpath(setting("dest"), logfile), "a") do f
                println(f, msg)
            end
        end
    end
    if category == 'i'
        logprint(msg)
    elseif category == 'd'
        setting("debug") && logprint("DEBUG: "*string(msg))
    elseif category == 'w'
        logprint("WARNING: "*string(msg), true)
    elseif category == 'e'
        logprint("ERROR: "*string(msg), true)
        exit(1)
    else
        simlog("Invalid log category "*string(category)*".", 'w')
    end
end

"""
    @simlog(message, category, args)

A wrapper macro around the function `simlog()`. Use this for debug statements, as it
prevents evaluation of the message string if debug mode is not on (a slight performance gain).
"""
macro simlog(msg::Union{String,Expr}, category::Char='i', args...)
    return :((($category != 'd') || setting("debug")) && simlog(:($$(esc(msg))), $category, $(args...)))
end

#XXX perhaps the utility functions below should go elsewhere?

"""
    diversity(w)

Calculate (average) alpha, beta and gamma diversity of the world.
Returns a tuple with the three values (a,b,g).
cf. Veech et al. 2002
"""
function diversity(world::Array{Patch,1})
    # calculate diversity with the Shannon-Wiener index
    function shannon(index::Dict{String,Int})
        isempty(index) && return 0
        total = sum(x -> index[x], keys(index))
        -1 * sum(s -> (index[s]/total)*log(index[s]/total), keys(index))
    end
    alphas = Float64[]
    globalindex = Dict{String, Int}()
    for p in world
        localindex = Dict{String, Int}()
        for s in unique(map(x -> x.lineage, p.community))
            localindex[s] = length(findall(x -> x.lineage == s, p.community))
        end
        merge!(+, globalindex, localindex)
        append!(alphas, shannon(localindex))
    end
    alpha = mean(alphas)
    gamma = shannon(globalindex)
    beta = gamma - alpha
    (alpha, beta, gamma)
end

"""
    heterozygosity(population)

Calculate the percentage of extraspecific chromosomes in a population's gene pool.
"""
function heterozygosity(population::Array{Individual,1})
    (!(setting("heterozygosity") || isempty(population))) && return 0
    lineage = population[1].lineage
    chromosomes = 0
    extraspecifics = 0
    for ind in population
        for c in ind.genome
            chromosomes += 1
            (c.lineage != lineage) && (extraspecifics += 1)
        end
    end
    return (extraspecifics/chromosomes) * 100
end

"""
    freespace(world)

Calculate the average amount of free space on each patch.
"""
function freespace(world::Array{Patch,1})
    space = 0
    for p in world
        if length(p.community) != 0
            space += p.capacity - sum(x -> x.size, p.community)
        else
            space += p.capacity
        end
    end
    round((space/length(world))/1e6, digits = 3)
end

"""
    idstring(individual)

A small utility function that returns a string identifier for a given individual.
"""
function idstring(ind::Individual)
    return ind.lineage*" "*string(ind.id)
end

max_stars_repo_name: CCTB-Ecomods/gemm

max_issues_repo_name: CCTB-Ecomods/gemm

max_forks_repo_name: CCTB-Ecomods/gemm


text:
using BioMedQuery.Entrez
using MySQL
using BioMedQuery.DBUtils

dbname="test"

config = Dict(:host=>"localhost", :dbname=>dbname, :username=>"root",
:pswd=>"", :overwrite=>true)

con = Entrez.DB.init_pubmed_db_mysql(config)

#check collection of tables
tables_query = BioMedQuery.DBUtils.select_all_tables(con)
tables = ["article","author","author2article","mesh_descriptor","mesh_heading","mesh_qualifier"]
@test sort(tables) == sort(tables_query)

#check minimum insert
BioMedQuery.DBUtils.insert_row!(con, "article", Dict(:pmid => 1234,
:title=>"Test Article",
:pubYear=>nothing))

#check duplicate error insert
duplicate_id = BioMedQuery.DBUtils.insert_row!(con, "article", Dict(:pmid => 1234,
    :title=>"Test Article", :pubYear=>nothing))

@test duplicate_id < 1
sel = BioMedQuery.DBUtils.db_select(con, ["pmid"], "article", Dict(:title=>"Test Article"))

@test length(sel[1]) == 1
@test sel[1][1] == 1234

#clean up
db_query(con, "DROP DATABASE IF EXISTS $dbname;")