text:
# g=5.0
# ENV["JULIA_CUDA_MEMORY_LIMIT"] = convert(Int, round(g * 1024 * 1024 * 1024))

include("exp.jl")

function main_ngraph()
    "Test how many graphs are needed to learn the model."
    # this requires the data generating spec to be different

    # TODO these are a lot of configs
    for d in [10, 20],
        ng in [100, 200, 500, 1000, 2000, 5000, 10000]

        specs = []
        for gtype in [:ER, :SF],
            k in [1]
            push!(specs, DataSpec(d=d, k=k, gtype=gtype,
                    noise=:Gaussian, mat=:CH3,
                    # when ng=100, total graph is 200, testing is 200*0.2=40, not enough for a batch of 64
                    ng=ng, bsize=32,
                    N=1))
        end
        specs = Array{DataSpec}(specs)

        @info "training .." d ng
        expID = exp_train(specs, eq2_model_fn,
                          # TODO I'll need to increase the training steps here
                          # CAUTION feed in the gtype in the model prefix
                          prefix="ngraph-ng=$ng-d=$d", train_steps=1e5,
                          test_throttle = 10,
                          merge=true)
    end
end

function main_ersf124()
    # test on different number of nodes
    for d in [10,20],
        (prefix, model_fn,nsteps) in [
            ("EQ2", eq2_model_fn, 1e5),
        ]

        specs = []
        for gtype in [:ER, :SF, :ER2, :SF2, :ER4, :SF4],
            k in [1]
            push!(specs, DataSpec(d=d, k=k, gtype=gtype,
                    noise=:Gaussian, mat=:CH3))
        end
        specs = Array{DataSpec}(specs)

        @info "training .." prefix d
        expID = exp_train(specs, model_fn,
                          # TODO I'll need to increase the training steps here
                          # CAUTION feed in the gtype in the model prefix
                          prefix="$prefix-ERSF124-d=$d", train_steps=nsteps,
                          test_throttle = 10,
                          merge=true)
    end
end

function main_data()
    # generate data for baseline usage
    for d in [10, 20, 50, 100, 150, 200, 300, 400]
        @info "Generating for" d ".."
        specs = []
        # this is combined specs. What was the filename?
        # will generate separate data files, and combine when loading
        for gtype in [:ER, :SF],
            k in [1]
            push!(specs, DataSpec(d=d, k=k, gtype=gtype,
                    noise=:Gaussian, mat=:CH3))
        end
        specs = Array{DataSpec}(specs)
        ds, test_ds = spec2ds(specs)
    end
end

function main_large()
    # large d, large ng and N, longer training time
    for (d, ng, N, bsize) in [
            (200, 1000, 1, 2),
            (150, 3000, 3, 4),
            # DEBUG test ng=10000
            (100, 10000, 3, 8),
            # DEBUG testing whether d=200 is trainable
            # d=200 is very slow
            # d=400
            (400, 1000, 1, 1)
        ],
        (prefix, model_fn,nsteps) in [
            # previous try didn't actually set ng, N, bsize, and it was marked
            # with old- prefix
            ("EQ2", eq2_model_fn, 5e4)
        ]

        specs = []
        for gtype in [:ER, :SF],
            k in [1]
            push!(specs, DataSpec(d=d, k=k, gtype=gtype,
                    noise=:Gaussian, mat=:CH3,
                    ng=ng, N=N, bsize=bsize))
        end
        specs = Array{DataSpec}(specs)

        # print more frequently for CNN and FC to get more data to print
        test_throttle = if prefix == "EQ2" 10 else 1 end

        @info "training .." prefix d
        expID = exp_train(specs, model_fn,
                          # TODO I'll need to increase the training steps here
                          # CAUTION feed in the gtype in the model prefix
                          prefix="$prefix-ERSF-k1-d=$d-ng=$ng-N=$N", train_steps=nsteps,
                          test_throttle = test_throttle,
                          merge=true)
    end
end

# Train on ER and test on SF and vice versa
function main_ersf()
    for d in [10,20,50, 100],
        (prefix, model_fn,nsteps) in [
            ("EQ2", eq2_model_fn, 3e4),
        ],
        gtype in [:ER, :SF]

        specs = DataSpec(d=d, k=1, gtype=gtype,
                        noise=:Gaussian, mat=:CH3)

        # print more frequently for CNN and FC to get more data to print
        test_throttle = if prefix == "EQ2" 10 else 1 end

        @info "training $prefix-$gtype-k1-d=$d .."
        expID = exp_train(specs, model_fn,
                          # TODO I'll need to increase the training steps here
                          # CAUTION feed in the gtype in the model prefix
                          prefix="$prefix-$gtype-k1-d=$d", train_steps=nsteps,
                          test_throttle = test_throttle,
                          merge=true)
        # Do the test here
        begin
            for gtype in [:ER, :SF]
                test_spec = DataSpec(d=d, k=1, gtype=gtype,
                                    noise=:Gaussian, mat=:CH3)
                exp_test(expID, test_spec, "TEST-gtype=$gtype")
            end
        end
    end
end

function main_ch3()
    for d in [10,20,50,100],
        (prefix, model_fn,nsteps) in [
            ("EQ2", eq2_model_fn, 5e4),
            ("FC", ()->fc_model(d=d, ch=2, z=1024, nlayer=6), 3e4),
            # ("FCreg", ()->fc_model(d=d, ch=2, z=1024, nlayer=6, reg=true), 3e4),
            ("CNN", ()->cnn_model(2), 3e4),
            # ("CNN2", ()->cnn_model(2, 128, (5,5), (2,2)), 3e4)
        ]

        specs = []
        for gtype in [:ER, :SF],
            k in [1]
            push!(specs, DataSpec(d=d, k=k, gtype=gtype,
                    noise=:Gaussian, mat=:CH3))
        end
        specs = Array{DataSpec}(specs)

        # print more frequently for CNN and FC to get more data to print
        test_throttle = if prefix == "EQ2" 10 else 1 end

        @info "training .." prefix d
        expID = exp_train(specs, model_fn,
                          # TODO I'll need to increase the training steps here
                          # CAUTION feed in the gtype in the model prefix
                          prefix="$prefix-ERSF-k1-d=$d", train_steps=nsteps,
                          test_throttle = test_throttle,
                          merge=true)
    end
end

function main_mat()
    # test the different matrix, including: CH3, COV, COR, ensemble on different k
    for d in [20,50],
        (mat, model_fn) in [
            (:CH3, eq2_model_fn),
            (:COV, eq_model_fn),
            (:COR, eq_model_fn)]

        specs = []
        for gtype in [:ER, :SF],
            k in [1,2,4,8]
            push!(specs, DataSpec(d=d, k=k, gtype=gtype,
                    noise=:Gaussian, mat=mat))
        end
        specs = Array{DataSpec}(specs)

        @info "training .." d mat
        expID = exp_train(specs, model_fn,
                          # TODO I'll need to increase the training steps here
                          # CAUTION feed in the gtype in the model prefix
                          prefix="EQ-d=$d-mat=$mat", train_steps=3e4,
                          merge=true)
    end
end

function main_ensemble_d()
    # ensemble on different d
    specs = []
    for gtype in [:ER, :SF],
        d in [20,30,40] # TODO use more aggressive d, e.g. 10, 30, 50

        push!(specs, DataSpec(d=d, k=1, gtype=gtype,
                noise=:Gaussian, mat=:CH3))
    end
    specs = Array{DataSpec}(specs)

    @info "training .."
    expID = exp_train(specs, eq2_model_fn,
                      # TODO I'll need to increase the training steps here
                      # CAUTION feed in the gtype in the model prefix
                      prefix="EQ2-CH3-d=[20,30,40]", train_steps=3e4,
                      merge=false)
end

function main_ensemble_d_cnn()
    # ensemble on different d
    specs = []
    for gtype in [:ER, :SF],
        d in [10,15,20] # TODO use more aggressive d, e.g. 10, 30, 50

        push!(specs, DataSpec(d=d, k=1, gtype=gtype,
                noise=:Gaussian, mat=:CH3))
    end
    specs = Array{DataSpec}(specs)

    @info "training .."
    expID = exp_train(specs, ()->cnn_model(2),
                      # TODO I'll need to increase the training steps here
                      # CAUTION feed in the gtype in the model prefix
                      prefix="CNN-CH3-d=[10,15,20]", train_steps=3e4,
                      test_throttle=1,
                      merge=false)
end

max_stars_repo_name: lihebi/DAG-EQ

max_issues_repo_name: lihebi/DAG-EQ

max_forks_repo_name: lihebi/DAG-EQ


text:
@testset "1610.maximum-number-of-visible-points.jl" begin
    @test max_visible_points([[2, 1], [2, 2], [3, 3]], 90, [1, 1]) == 3
    @test max_visible_points([[2, 1], [2, 2], [3, 4], [1, 1]], 90, [1, 1]) == 4
    @test max_visible_points([[1, 0], [2, 1]], 13, [1, 1]) == 1
end

ext: jl


text:
module LDA_SIM

using FileIO, JLD2, Turing, ContinuousBenchmarks

# V <- 5; # words: river, stream, bank, money, loan
# K <- 2; # topics: RIVER, BANK
# M <- 25;  # docs
V = 100
K = 5
M = 10

# alpha <- rep(1/K,K);
# beta <- rep(1/V,V);
alpha = collect(ones(K) / K);
beta = collect(ones(V) / V);

# phi <- array(NA,c(2,5));
# phi[1,] = c(0.330, 0.330, 0.330, 0.005, 0.005);
# phi[2,] = c(0.005, 0.005, 0.330, 0.330, 0.330);
phi = rand(Dirichlet(beta), K)
#println(phi)

# theta <- rdirichlet(M,alpha);
theta = rand(Dirichlet(alpha), M)
#println(theta)

# avg_doc_length <- 10;
avg_doc_length = 1000

# doc_length <- rpois(M,avg_doc_length);
doc_length = rand(Poisson(avg_doc_length), M)

# N <- sum(doc_length);
N = sum(doc_length)

# w <- rep(NA,N);
# doc <- rep(NA,N);
w = Vector{Int}(undef, N)
doc = Vector{Int}(undef, N)

# n <- 1;
# for (m in 1:M) {
#   for (i in 1:doc_length[m]) {
#     z <- which(rmultinom(1,1,theta[m,]) == 1);
#     w[n] <- which(rmultinom(1,1,phi[z,]) == 1);
#     doc[n] <- m;
#     n <- n + 1;
#   }
# }
n = 1
for m = 1:M, i=1:doc_length[m]
    global n
    z = rand(Categorical(theta[:,m]))
    w[n] = rand(Categorical(phi[:,z]))
    doc[n] = m
    n = n + 1
end

ldastandata = [
    Dict(
    "K" => K,
    "V" => V,
    "M" => M,
    "N" => N,
    "w" => w,
    "doc" => doc,
    "alpha" => alpha,
    "beta" => beta,
    # "phi" => phi,
    # "theta" => theta,
    # "D" => M * K + K * V
  )
]

# println(ldastandata)

save(joinpath(ContinuousBenchmarks.STAN_DATA_DIR, "ldastandata.jld2"), "data", ldastandata)

end


text:
using Test
using LogicCircuits
using ProbabilisticCircuits
using LightGraphs
using MetaGraphs
using Suppressor

@testset "Load a small PSDD and test methods" begin
   file = zoo_psdd_file("little_4var.psdd")
   prob_circuit = load_prob_circuit(file);
   @test prob_circuit isa ProbCircuit

   # Testing number of nodes and parameters
   @test  9 == num_parameters(prob_circuit)
   @test 20 == num_nodes(prob_circuit)

   # Testing Read Parameters
   EPS = 1e-7
   or1 = children(children(prob_circuit)[1])[2]
   @test abs(or1.log_probs[1] - (-1.6094379124341003)) < EPS
   @test abs(or1.log_probs[2] - (-1.2039728043259361)) < EPS
   @test abs(or1.log_probs[3] - (-0.916290731874155)) < EPS
   @test abs(or1.log_probs[4] - (-2.3025850929940455)) < EPS

   or2 = children(children(prob_circuit)[1])[1]
   @test abs(or2.log_probs[1] - (-2.3025850929940455)) < EPS
   @test abs(or2.log_probs[2] - (-2.3025850929940455)) < EPS
   @test abs(or2.log_probs[3] - (-2.3025850929940455)) < EPS
   @test abs(or2.log_probs[4] - (-0.35667494393873245)) < EPS

   @test abs(prob_circuit.log_probs[1] - (0.0)) < EPS
end

psdd_files = ["little_4var.psdd", "msnbc-yitao-a.psdd", "msnbc-yitao-b.psdd", "msnbc-yitao-c.psdd", "msnbc-yitao-d.psdd", "msnbc-yitao-e.psdd", "mnist-antonio.psdd"]

@testset "Test parameter integrity of loaded PSDDs" begin
   for psdd_file in psdd_files
      @test check_parameter_integrity(load_prob_circuit(zoo_psdd_file(psdd_file)))
   end
end

@testset "Test parameter integrity of loaded structured PSDDs" begin
   circuit, vtree = load_struct_prob_circuit(
      zoo_psdd_file("little_4var.psdd"), zoo_vtree_file("little_4var.vtree"))
   @test check_parameter_integrity(circuit)
   @test vtree isa PlainVtree
   @test respects_vtree(circuit, vtree)
end

@testset "Test loaded logistic circuits" begin
   classes = 2
   lc = zoo_lc("little_4var.circuit", classes)
   @test lc isa LogisticCircuit
   @test num_nodes(lc) == 29
   @test num_parameters(lc) == num_parameters_per_class(lc) * classes
   @test check_parameter_integrity(lc)
   @test all(lc.thetas ≈ [0.6803363307333976 0.7979990493834855])
   @test all(lc.children[1].thetas ≈ [0.34978385969447856 0.8620937209951014])
   or1 = lc.children[1].children[1].children[1]
   @test all(or1.thetas ≈ [0.012723368087276588 0.44600859247854274;
                           0.4126851950485019 0.8204476705654302;
                           0.2633016148946008 0.010045227037839832;
                           0.27132714432580674 0.9677387544772587])
end

@testset "Test loaded clt file" begin
   # TODO add more clt files
   t = zoo_clt("4.clt")
   @test t isa CLT
   @test_nowarn @suppress_out print_tree(t)
   @test nv(t) == 4
   @test ne(t) == 0
   for (v, p) in enumerate([.2, .1, .3, .4])
      cpt = get_prop(t, v, :cpt)
      @test cpt[1] + cpt[0] ≈ 1.0
      @test cpt[1] ≈ p
   end
end


text:
@test mdinfo() !== nothing
groups = ["symmetric", "inverse", "illcond", "posdef", "eigen","sparse", "random", "regprob", "all"]

for group in groups
    @test !isempty(listnames(Symbol(group)))
end

@test_throws DataError matrixdepot("something")

n = rand(1:55)
name = mdlist(builtin(n))
@test mdinfo(name) !== nothing

nlist = mdlist(builtin([1, 3]) | builtin(4:20))
m = length(mdlist("**"))
@test isempty(mdinfo(builtin(m+1)))

@addgroup newlist = mdlist(builtin(3:6) | builtin(20))

@test mdlist(:newlist) == mdlist(builtin(3:6,20))

@rmgroup newlist

# testing the new API
#
# it is assumed that all matrices are generated during "test/download.jl"
# and the builtin- and user-defined "randsym" but no others are available
#

@testset "mdlist" begin

REM = length(mdlist("*/*"))
@test REM >= 2500
#@test REM in [2757, 2856]   # depends on whether UFL or TAMU url has been used
@test length(mdlist(:builtin)) == 59
@test length(mdlist(:user)) in [0, 1]

@test mdlist("") == []
@test mdlist("HB/1138_bus") == ["HB/1138_bus"]
@test mdlist(sp(1)) == ["HB/1138_bus"]
@test mdlist(mm(1)) == ["Harwell-Boeing/psadmit/1138_bus"]
@test sort(mdlist(sp(1:3000))) == mdlist("*/*")
@test sort(mdlist(mm(1:3000))) == mdlist("*/*/*")
@test mdlist(builtin(:)) == mdlist(isbuiltin)
@test mdlist(user(:)) == mdlist(isuser)
@test mdlist("*") == mdlist(islocal)
@test length(mdlist("HB/*")) == 292
@test mdlist("HB**") == mdlist("HB/**")
@test_throws ArgumentError  listdir("*/*")
@test listdir("Harwell-Boeing//") == ["Harwell-Boeing/*/* - (292)"]
@test mdlist("*/hamm/*") == ["misc/hamm/add20", "misc/hamm/add32", "misc/hamm/memplus"]
@test mdlist("*/hamm/*a*3?") == ["misc/hamm/add32"]
@test mdlist("*/hamm/*a*3[123]") == ["misc/hamm/add32"]
@test mdlist("*/hamm/*a*3[123]?") == []
@test length(mdlist("*/*/*")) == 498

@test listdir("*//*/*") == ["Harwell-Boeing/*/* - (292)", "NEP/*/* - (73)",
                           "SPARSKIT/*/* - (107)", "misc/*/* - (26)"]
@test listdir("//*") == ["/* - ($(length(mdlist(:local))))"]
@test listdir("//*/*") == ["/*/* - ($REM)"]
@test listdir("//*/*/*") == ["/*/*/* - (498)"]
@test listdir("HB/") == ["HB/* - (292)"]
@test length(mdlist("Harwell-Boeing/*/*")) == 292
@test mdlist(r".*ng/ma.*") == ["Harwell-Boeing/manteuffel/man_5976"]
@test mdlist(sp(2001:2002)) == ["JGD_Groebner/c8_mat11_I", "JGD_Groebner/f855_mat9"]
@test length(mdlist(sp(2757:3000))) == REM - 2756
@test_throws ArgumentError mdlist(:xxx)
@test length(mdlist(isremote)) == REM + 498
@test length(mdlist(isloaded)) + length(mdlist(isunloaded)) == length(mdlist(isremote))
@test length(mdlist(isbuiltin)) + length(mdlist(isuser)) == length(mdlist(islocal))
@test length(mdlist(islocal)) + length(mdlist(isremote)) == length(mdlist(:all))
@test length(mdlist("**")) == length(mdlist("*")) + length(mdlist("*/*")) + length(mdlist("*/*/*"))
@test mdlist(:all) == mdlist("**")
@test length(mdlist(:symmetric)) == 22
@test length(mdlist(:illcond)) == 20

# intersections and unions
@test mdlist((:posdef, :sparse)) == ["poisson", "wathen"]
@test length(mdlist([:posdef,:sparse])) + length(mdlist((:posdef,:sparse))) == length(mdlist(:posdef)) + length(mdlist(:sparse))

# predicates of remote and local matrices
@test length(mdlist(issymmetric)) == 30
@test length(mdlist(:symmetric & "kahan")) == 0
@test length(mdlist(:symmetric & "hankel")) == 1
@test mdlist(:local) == mdlist(islocal)
@test mdlist(:builtin) == mdlist(isbuiltin)
@test mdlist(:user) == mdlist(isuser)
@test mdlist(:symmetric) == mdlist(issymmetric & islocal)

# general metadata
@test length(mdlist(@pred(m < 10000) & "HB/*")) == 284   # items with m < *
@test length(mdlist(@pred(m < 10000) & isloaded)) == 10   # items with m < *
@test length(mdlist(@pred(n < 10000) & "HB/*")) == 283   # items with n < *
@test length(mdlist(@pred(n < 10000) &  isloaded)) == 9    # items with n < *
@test length(mdlist(@pred(nnz < 5000) & "HB/*")) == 163    # items with nnz < *
@test length(mdlist(@pred(nnz < 5000) & isloaded)) == 2    # items with nnz < *
@test length(mdlist(@pred(m > n) & "HB/*")) == 9   # items with m > n
@test length(mdlist(@pred(m > n) & isloaded)) == 0   # items with m > n

# metadata from header files
@test length(mdlist(@pred(occursin("problem", kind)) & isloaded)) == 6
@test length(mdlist(keyword("graph") & isloaded)) == 2
@test length(mdlist(@pred(0 < date <= 1990) & isloaded)) == 1
@test length(mdlist(@pred(date >= 2006) & @pred(0 < date <= 2006) & isloaded)) == 2
@test length(mdlist(@pred(date == 0) & mm(:) & isloaded)) == 3

# metadata from ss_index.mat
@test length(mdlist(@pred(posdef) & "HB/*")) == 60
@test length(mdlist(isposdef)) >= 60
@test mdlist(@pred(posdef)) == mdlist(isposdef & isremote)
@test mdlist(:posdef) == mdlist(isposdef & islocal)

# metadata from svd files
@test mdlist(@pred(svdok)) == ["HB/1138_bus"]
data = mdopen("HB/1138_bus").data
@test length(data.sv) > 0
@test data.cond > 8.0e6
@test data.norm > 30000
@test data.rank == 1138
@test data.svgap == Inf
@test data.cholcand

end

@testset "logical" begin
# for the boolean syntax
@test mdlist(~islocal) == mdlist(isremote)
@test length(mdlist(isloaded & issymmetric)) == 7
@test length(mdlist(isloaded & ~issymmetric)) == 4
@test length(mdlist(isloaded & ~issymmetric | isuser & issymmetric)) == 5
@test length(mdlist(!islocal & issymmetric | isuser & issymmetric)) == 9
@test mdlist(islocal & ~isbuiltin) == mdlist(isuser)
@test mdlist(islocal & ~isbuiltin) == mdlist(isuser)
@test "a" & "b" === ("a", "b")
@test ~"a" & "b" === (~"a", "b")
@test ~"a" * "b" === ~"ab"
@test ~"ab" === ~'a' * 'b'
@test ~~islocal === islocal
@test mdlist(~"*a*" & ~ "*e*") == mdlist(~["*a*", "*e*"])
@test mdlist(()) == mdlist(:all)
@test "a" & r"b" == ("a", r"b")
@test "a" & r"b" & "c" == ("a", r"b", "c")
@test "a" | r"b" == ["a", r"b"]
@test "a" | r"b" | "c" == ["a", r"b", "c"]
@test mdlist(islocal & ~("*a*" | "*e*")) == mdlist(:local & ~"*a*" & ~ "*e*")
@test "a" & ( "b" & "c" ) == ("a", "b", "c")
@test "a" | ( "b" | "c" ) == ["a", "b", "c"]
@test_throws ArgumentError mdlist([] & :invalid_group_name)

@test MatrixDepot.mdlist(builtin(10, 1:7, 3)) == ["baart", "binomial", "blur", "cauchy",
                                                "chebspec", "chow", "circul", "deriv2"]

@test MatrixDepot.fname(sin) == "unknown-function"
@test_throws MethodError mdopen("baart", 10, 11)
end

@testset "properties" begin

md = mdopen("gravity", 10)
io = IOBuffer()
@test show(io, md) === nothing
@test propertynames(md) == [:A, :m, :n, :nnz, :dnz]
@test md.A isa AbstractMatrix
@test md.nnz == count(md.A .!= 0)
@test md.dnz == 0
@test_throws DataError md.x

md = mdopen("gravity", 10, false)
io = IOBuffer()
@test show(io, md) === nothing
@test propertynames(md) == [:A, :b, :x, :m, :n, :nnz, :dnz]
@test md.A isa AbstractMatrix
@test md.b isa AbstractVector
@test md.x isa AbstractVector

md = mdopen("HB/well1033")
io = IOBuffer()
@test show(io, md) === nothing
@test propertynames(md) == [:A, :b, :m, :n, :nnz, :dnz]
@test md.A isa AbstractMatrix
@test md.b isa AbstractMatrix
@test_throws DataError md.x
@test_throws DataError md.invalidname
@test_throws ErrorException md.data.invalidname
@test md.dnz == md.data.header.nnz

@test length(mdlist(keyword("graph" & ("butterfly" | "network")) & isloaded)) == 2
@test length(mdlist(hasdata(:A & (:b | :x)))) == 2

end


text:
#*******************************************************/
#* Copyright(c) 2018 by Artelys                        */
#* This source code is subject to the terms of the     */
#* MIT Expat License (see LICENSE.md)                  */
#*******************************************************/
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This file contains routines to implement problemDef.h for
# a simple nonlinear least-squares problem.
#
# min  ( x0*1.309^x1 - 2.138 )^2 +( x0*1.471^x1 - 3.421 )^2 +( x0*1.49^x1 - 3.597 )^2
#        +( x0*1.565^x1 - 4.34 )^2 +( x0*1.611^x1 - 4.882 )^2 +( x0*1.68^x1-5.66 )^2
#
# The standard start point(1.0, 5.0) usually converges to the standard
# minimum at(0.76886, 3.86041), with final objective = 0.00216.
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

using KNITRO
using Test

function example_lsq2(; verbose=true)
    #*------------------------------------------------------------------*
    #*     FUNCTION callbackEvalR                                       *
    #*------------------------------------------------------------------*
    # The signature of this function matches KNITRO.KN_eval_callback in knitro.h.
    # Only "rsd" is set in the KNITRO.KN_eval_result structure.
    function callbackEvalR(kc, cb, evalRequest, evalResult, userParams)
        if evalRequest.evalRequestCode != KNITRO.KN_RC_EVALR
            println(
                "*** callbackEvalR incorrectly called with eval type ",
                evalRequest.evalRequestCode,
            )
            return -1
        end
        x = evalRequest.x

        # Clamp x1 to 1000 so as to avoid overflow
        if x[2] > 1000.0
            x[2] = 1000.0
        end

        # Evaluate nonlinear residual components
        evalResult.rsd[1] = x[1] * 1.309^x[2]
        evalResult.rsd[2] = x[1] * 1.471^x[2]
        evalResult.rsd[3] = x[1] * 1.49^x[2]
        evalResult.rsd[4] = x[1] * 1.565^x[2]
        evalResult.rsd[5] = x[1] * 1.611^x[2]
        evalResult.rsd[6] = x[1] * 1.68^x[2]

        return 0
    end

    #*------------------------------------------------------------------*
    #*     FUNCTION callbackEvalRJ                                      *
    #*------------------------------------------------------------------*
    # The signature of this function matches KNITRO.KN_eval_callback in knitro.h.
    # Only "rsdJac" is set in the KNITRO.KN_eval_result structure.
    function callbackEvalRJ(kc, cb, evalRequest, evalResult, userParams)
        if evalRequest.evalRequestCode != KNITRO.KN_RC_EVALRJ
            println(
                "*** callbackEvalRJ incorrectly called with eval type %d" %
                evalRequest.evalRequestCode,
            )
            return -1
        end
        x = evalRequest.x

        # Clamp x1 to 1000 so as to avoid overflow
        if x[2] > 1000.0
            x[2] = 1000.0
        end

        # Evaluate non-zero residual Jacobian elements(row major order).
        evalResult.rsdJac[1] = 1.309^x[2]
        evalResult.rsdJac[2] = x[1] * log(1.309) * 1.309^x[2]
        evalResult.rsdJac[3] = 1.471^x[2]
        evalResult.rsdJac[4] = x[1] * log(1.471) * 1.471^x[2]
        evalResult.rsdJac[5] = 1.49^x[2]
        evalResult.rsdJac[6] = x[1] * log(1.49) * 1.49^x[2]
        evalResult.rsdJac[7] = 1.565^x[2]
        evalResult.rsdJac[8] = x[1] * log(1.565) * 1.565^x[2]
        evalResult.rsdJac[9] = 1.611^x[2]
        evalResult.rsdJac[10] = x[1] * log(1.611) * 1.611^x[2]
        evalResult.rsdJac[11] = 1.68^x[2]
        evalResult.rsdJac[12] = x[1] * log(1.68) * 1.68^x[2]

        return 0
    end

    #*------------------------------------------------------------------*
    #*     main                                                         *
    #*------------------------------------------------------------------*

    # Create a new Knitro solver instance.
    kc = KNITRO.KN_new()

    # Add the variables/parameters.
    # Note: Any unset lower bounds are assumed to be
    # unbounded below and any unset upper bounds are
    # assumed to be unbounded above.
    n = 2 # # of variables/parameters
    KNITRO.KN_add_vars(kc, n)

    # In order to prevent the possiblity of numerical
    # overflow from very large numbers, we set a
    # reasonable upper bound on variable x[1] and set the
    # "honorbnds" option for this variable to enforce
    # that all trial x[1] values satisfy this bound.
    KNITRO.KN_set_var_upbnd(kc, 1, 100.0)
    KNITRO.KN_set_var_honorbnd(kc, 1, KNITRO.KN_HONORBNDS_ALWAYS)

    # Add the residuals.
    m = 6 # # of residuals
    KNITRO.KN_add_rsds(kc, m)

    # Set the array of constants in the residuals
    KNITRO.KN_add_rsd_constants_all(kc, [-2.138, -3.421, -3.597, -4.34, -4.882, -5.66])

    # Add a callback function "callbackEvalR" to evaluate the nonlinear
    # residual components.  Note that the constant terms are added
    # separately above, and will not be included in the callback.
    cb = KNITRO.KN_add_lsq_eval_callback(kc, callbackEvalR)

    # Also add a callback function "callbackEvalRJ" to evaluate the
    # Jacobian of the residuals.  If not provided, Knitro will approximate
    # the residual Jacobian using finite-differencing.  However, we recommend
    # providing callbacks to evaluate the exact Jacobian whenever
    # possible as this can drastically improve the performance of Knitro.
    # We specify the residual Jacobian in "dense" row major form for simplicity.
    # However for models with many sparse residuals, it is important to specify
    # the non-zero sparsity structure of the residual Jacobian for efficiency
    #(this is true even when using finite-difference gradients).
    KNITRO.KN_set_cb_rsd_jac(kc, cb, KNITRO.KN_DENSE_ROWMAJOR, callbackEvalRJ)

    kn_outlev = verbose ? KNITRO.KN_OUTLEV_ALL : KNITRO.KN_OUTLEV_NONE
    KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_OUTLEV, kn_outlev)

    # Solve the problem.
    #
    # Return status codes are defined in "knitro.h" and described
    # in the Knitro manual.

    nRC = KNITRO.KN_solve(kc)

    # An example of obtaining solution information.
    # Return status codes are defined in "knitro.h" and described
    # in the Knitro manual.
    nStatus, obj, x, lambda_ = KNITRO.KN_get_solution(kc)

    if nStatus != 0 && verbose
        println("Knitro successful. The optimal solution is:")
        for i in 1:n
            println("x[$i]= ", x[i])
        end
    end

    # Delete the knitro solver instance.
    KNITRO.KN_free(kc)

    @testset "Example LSQ2" begin
        @test nStatus == 0
        @test obj ≈ 0.00216 atol = 1e-4
        @test x ≈ [0.76886, 3.86041] atol = 1e-4
    end
end

example_lsq2(; verbose=isdefined(Main, :KN_VERBOSE) ? KN_VERBOSE : true)


text:
################################################################################
# COMPLEX BLOCKS
################################################################################

"""
    build_processor_tile(num_links, kwargs)

Build a processor tile.
"""
function build_processor_tile(
        style;
        include_memory = false,
        name = include_memory ? "memory_processor_tile" : "processor_tile",
    )

    num_fifos = 2
    num_links = links(style)
    # No need to assign metadata to this component.
    component = Component(name)

    # General metadata for ciruit link routing components.
    cl_metadata = routing_metadata("circuit_link")

    # Add the circuit switched ports
    for dir in directions(style)
        for (suffix,class)  in zip(("_in", "_out"), (Input, Output))
            meta = top_level_port_metadata(style, dir, class, "circuit_link", num_links)
            port_name = join((dir, suffix))
            add_port(component, port_name, class, num_links, metadata = meta)
        end
    end

    # Instantiate the processor primitive - do special assignment for ASAP mapping.
    proc_component = build_processor(style, include_memory = include_memory)

    add_child(component, proc_component, "processor")

    # Instantiate the directional routing muxes
    routing_mux = build_mux(length(directions(style)),1)
    for dir in directions(style)
        name = "$(dir)_mux"
        add_child(component, routing_mux, name, num_links)
    end

    # Instantiate the muxes routing data to the fifos
    # Allocate a mux input for each input on the circuit net.
    num_fifo_inputs = length(directions(style)) * num_links
    fifo_mux = build_mux(num_fifo_inputs, 1, metadata = cl_metadata)
    add_child(component, fifo_mux, "fifo_mux", num_fifos)

    # Add memory ports and links if necessary.
    if include_memory
        return_link  = routing_metadata("memory_response_link")
        request_link = routing_metadata("memory_request_link")
        add_port(component, "memory_in", Input, metadata = return_link)
        add_port(component, "memory_out", Output, metadata = request_link)
        add_link(component, "processor.memory_out", "memory_out", metadata = request_link)
        add_link(component, "memory_in", "processor.memory_in", metadata = return_link)
    end

    # Connect outputs of muxes to the tile outputs
    for dir in directions(style), i = 0:num_links-1
        # Each mux has only one output.
        mux_port    = "$(dir)_mux[$i].out[0]"
        tile_port   = "$(dir)_out[$i]"
        add_link(component, mux_port, tile_port, metadata = cl_metadata)
    end

    # Circuit switch output links
    for dir in directions(style), i = 0:num_links-1
        # Make the name for the processor.
        proc_port = "processor.$dir[$i]"
        mux_port = "$(dir)_mux[$i].in[0]"
        add_link(component, proc_port, mux_port, metadata = cl_metadata)
    end

    # Connect input fifos.
    for i = 0:num_fifos-1
        fifo_port = "fifo_mux[$i].out[0]"
        proc_port = "processor.fifo[$i]"
        add_link(component, fifo_port, proc_port, metadata = cl_metadata)
    end

    # Connect input ports to inputs of muxes

    # Tracker Structure storing what ports on a mux have been used.
    # Initialize to 1 because the processor is connect to port 0 on all
    # multiplexors.
    index_tracker = Dict((d,i) => 1 for d in directions(style), i in 0:num_links)
    # Add entries for the fifo
    for i in 0:num_fifos
        index_tracker[("fifo",i)] = 0
    end
    for dir in directions(style)
        for i in 0:num_links-1
            # Create a source port for the tile input.
            source_port = "$(dir)_in[$i]"
            # Begin adding sink ports.
            # Need one for each mux on this layer and num_fifos for each of
            # the input fifos.
            sink_ports = String[]
            # Go through all fifos that directions that are not the current
            # directions.
            for d in Iterators.filter(!isequal(dir), directions(style))
                # Get the next free port for this mux.
                key = (d,i)
                index = index_tracker[key]
                mux_port = "$(d)_mux[$i].in[$index]"
                # Add this port to the list of sink ports
                push!(sink_ports, mux_port)
                # increment the index tracker
                index_tracker[key] += 1
            end
            # Add the fifo mux entries.
            for j = 0:num_fifos-1
                key = ("fifo", j)
                index = index_tracker[key]
                fifo_port = "fifo_mux[$j].in[$index]"
                push!(sink_ports, fifo_port)
                index_tracker[key] += 1
            end
            # Add the connection to the component.
            add_link(component, source_port, sink_ports, metadata = cl_metadata)
        end
    end
    check(component)
    return component
end

################################################################################
#                           processor
################################################################################
function build_processor(
        style;
        include_memory = false,
        num_fifos = 2,
    )
    num_links = links(style)
    # Build the metadata dictionary for the processor component
    metadata = Dict{String,Any}()

    if include_memory
        comp_metadata = mem_proc_metadata()
        name = "memory_processor"
    else
        comp_metadata = proc_metadata()
        name = "standard_processor"
    end

    component = Component(name, metadata = comp_metadata)
    # fifos
    add_port(component, "fifo", Input, num_fifos, metadata = proc_fifo_metadata(num_fifos))
    # ports. Neet to play some indexing games for the metadata vector to get
    # indices to line up with the Asap4 manual.
    port_metadata = proc_output_metadata(length(directions(style)), num_links)
    # NOTE: Iterators.product iterates over the first element the quickest.
    for (port_index,(str,i)) in enumerate(Iterators.product(directions(style), 1:num_links))
        add_port(component, "$str[$(i-1)]", Output, metadata = port_metadata[port_index])
    end
    # Add memory ports. Will only be connected in the memory processor tile.
    if include_memory
        add_port(component, "memory_in", Input, metadata = proc_memory_return_metadata())
        add_port(component, "memory_out", Output, metadata = proc_memory_request_metadata())
    end
    # Return the created type
    return component
end

function build_memory(nports = 2)
    component = Component("memory_$(nports)port", metadata = mem_nport_metadata(nports))

    # Instantiate ports
    add_port(component, "in", Input, nports, metadata = mem_memory_request_metadata(nports))
    add_port(component, "out", Output, nports, metadata = mem_memory_return_metadata(nports))

    return component
end

##############################
#       INPUT HANDLER        #
##############################
build_input_handler(style::Style) = build_input_handler(iolinks(style))
function build_input_handler(num_links::Integer)
    component = Component("input_handler", metadata = input_handler_metadata())
    add_port(
        component,
        "out",
        Output,
        num_links,
        metadata = input_handler_port_metadata(num_links)
    )

    return component
end

##############################
#       OUTPUT HANDLER       #
##############################
build_output_handler(style::Style) = build_output_handler(iolinks(style))
function build_output_handler(num_links::Integer)
    component = Component("output_handler", metadata = output_handler_metadata())
    add_port(
        component,
        "in",
        Input,
        num_links,
        metadata = output_handler_port_metadata(num_links)
    )

    return component
end

################################################################################
# Functions for connecting processors, IO, and memories together
################################################################################
squash(x) = reshape(x, :)
function connect_processors(toplevel, style)

    # Build metadata dictionary for capacity and cost
    metadata = Dict(
        "cost"          => 1.0,
        "capacity"      => 1,
        "link_class"    => "circuit_link"
    )

    for rule in procrules(style)
        connection_rule(toplevel, rule, metadata = metadata)
    end
    return nothing
end

function connect_io(toplevel, style)

    # Build metadata dictionary for capacity and cost
    metadata = Dict(
        "cost"          => 1.0,
        "capacity"      => 1,
        "link_class"    => "link_class"
    )

    for rule in iorules(style)
        connection_rule(toplevel, rule, metadata = metadata)
    end
    return nothing
end

function connect_memories(toplevel, style)
    # Create metadata dictionary for the memory links.
    request_metadata = Dict(
        "cost"          => 1.0,
        "capacity"      => 1,
        "link_class"    => "memory_request_link",
   )

    return_metadata = Dict(
        "cost"          => 1.0,
        "capacity"      => 1,
        "link_class"    => "memory_response_link",
   )
    ###########################
    # Connect 2 port memories #
    ###########################

    for rule in memory_request_rules(style)
        connection_rule(toplevel, rule, metadata = request_metadata)
    end

    for rule in memory_return_rules(style)
        connection_rule(toplevel, rule, metadata = return_metadata)
    end

    return nothing
end


text:
using MathProgBase, Gurobi
using Distributions, DataFrames
using Compat

InputVector{T<:Union{Real,Char}} = Union{Vector{T},Real,Char}
const SymbolInputVector = Union{Vector{Symbol},Symbol}

"""
Builds a mixed integer programming problem as defined in MathProgBase.
Args:
     c:      is the objective vector, always in the sense of minimization
     A:      is the constraint matrix
     sense:  is a vector of constraint sense characters '<', '=', and '>'
     b:      is the right-hand side vector
     l:      is the vector of lower bounds on the variables
     u:      is the vector of upper bounds on the variables, and
     solver: specifies the desired solver, see :ref:`choosing solvers <choosing-solvers>`.

A scalar is accepted for the b, sense, l, and u arguments, in which case its value is replicated.
The values -Inf and Inf are interpreted to mean that there is no corresponding lower or upper bound.
"""
function buildmio(c::InputVector, A::AbstractMatrix, rowlb::InputVector, rowub::InputVector, vartypes::SymbolInputVector, lb::InputVector, ub::InputVector, solver)
    m = MathProgBase.LinearQuadraticModel(solver)
    nrow,ncol = size(A)

    c = MathProgBase.HighLevelInterface.expandvec(c, ncol)
    rowlbtmp = MathProgBase.HighLevelInterface.expandvec(rowlb, nrow)
    rowubtmp = MathProgBase.HighLevelInterface.expandvec(rowub, nrow)
    lb = MathProgBase.HighLevelInterface.expandvec(lb, ncol)
    ub = MathProgBase.HighLevelInterface.expandvec(ub, ncol)
    vartypes = MathProgBase.HighLevelInterface.expandvec(vartypes, ncol)

    # rowlb is allowed to be vector of senses
    if eltype(rowlbtmp) == Char
        realtype = eltype(rowubtmp)
        sense = rowlbtmp
        rhs = rowubtmp
        @assert realtype <: Real
        rowlb = Array{realtype}(undef, nrow)
        rowub = Array{realtype}(undef, nrow)
        for i in 1:nrow
            if sense[i] == '<'
                rowlb[i] = typemin(realtype)
                rowub[i] = rhs[i]
            elseif sense[i] == '>'
                rowlb[i] = rhs[i]
                rowub[i] = typemax(realtype)
            elseif sense[i] == '='
                rowlb[i] = rhs[i]
                rowub[i] = rhs[i]
            else
                error("Unrecognized sense '$(sense[i])'")
            end
        end
    else
        rowlb = rowlbtmp
        rowub = rowubtmp
    end

    MathProgBase.loadproblem!(m, A, lb, ub, c, rowlb, rowub, :Min)
    MathProgBase.setvartype!(m, vartypes)
    return m
end

"""
Solves the maximization problem max |beta0*x(i,1)+x(i,:)*t| over t subject to bnd.
Args:
    x (array):      (n by k) matrix of covariate data.
    beta0 (int):    coefficient for the first covariate in x.
    bnd (array):    ((k-1) by 2) matrix where the 1st and 2nd columns store the lower and
                    upper bounds of the unknown coefficients.
"""
function MIObnd(x::AbstractArray,beta0::Integer,bnd::AbstractArray)
    println("Start mio bound")
    n=size(x)[1]
    k=size(x)[2]-1

    v = zeros(2,1)
    value=zeros(n,1)

    #Set Params
    tol=0.000001
    solver = GurobiSolver(OutputFlag=0,OptimalityTol=tol,FeasibilityTol=tol,IntFeasTol=tol)

    for i = 1:n
        println("i = ",i)
        alpha = beta0*x[i,1]
        obj = x[i,2:k+1]

        try
            A = x[i,2:k+1][:,:]'
            rhs = -alpha
            result= linprog(-obj, A, '>', rhs, bnd[:,1], bnd[:,2], solver)
            v[1]=-result.objval+alpha
            catch
                println("Error Reported")
        end

        obj = -x[i,2:k+1]
        try
            A = -x[i,2:k+1][:,:]'
            rhs = alpha
            result = linprog(-obj, A, '>', rhs, bnd[:,1], bnd[:,2], solver)
            v[2] = -result.objval-alpha
            catch
                println("Error Reported")
        end
        value[i]=maximum(v)
    end
    return value
end

"""
Calculates the maximum score of the objective function.
Args:
    y (array):      vector of binary outcomes.
    x_foc (array):  (n by k) matrix of focus covariates where the first column should
                    contain data of the continous regressor with respect to the scale normalization.
    x_aux (array):  (n by d) matrix of auxiliary covariates which will be selected based on best
                    subset covariate selection procedure.
    beta0 (int):    Coefficient taking value either 1 or -1 to normalize the scale for the first
                    first covariate in x_foc.
    q (int):        Cardinality constraint for the covariate selection.
    abgap (int):    Absolute gap specified for early termination of the MIO solver.
    bnd (array):    (((k-1)+d) by 2) matrix where the first and second columns respectively
                    store the lower and upper bounds of the unknown coefficients. First (k-1)
                    rows correspond to the bounds of the focused covariates excluding the first one.
                    Remaining d rows correspond to the bounds of the auxiliary covariates.
    mio (int) :     MIO formulation of the best subset maximum score problem
                    set mio = 1 for Method 1 of the paper
                    set mio = 2 for Method 2 of the paper
Yields:
    bhat (array):   Maximum score estimates for the unknown coefficients.
    score (array):  Value of maximum score objective function.
    gap   (array):  MIO optimization gap value in case of early termination
                    gap = 0 ==> optimal solution is found within the time limit.
    rtime (array):  Time used by the MIO solver in the estimation procedure.
"""
function MaxScore(y::AbstractArray,x_foc::AbstractArray,x_aux::AbstractArray,beta0::Integer,q::Integer,T::Real,abgap::Real,bnd::AbstractArray,mio::Integer)
    N=length(y)
    k=size(x_foc)[2]-1
    d=size(x_aux)[2]
    bhat=zeros(k+d)
    score = 0
    gap = 0
    rtime = 0
    ncount = 0

    miobnd=MIObnd([x_foc x_aux],beta0,bnd)

    #Set Params
    tol=0.000001
    env = Gurobi.Env()
    setparams!(env; OutputFlag=0,OptimalityTol=tol,FeasibilityTol=tol,IntFeasTol=tol)
    solver = GurobiSolver(env)

    if T > 0
    setparams!(env; TimeLimit = T)
    solver = GurobiSolver(env)
    end

    if abgap > 0
    setparams!(env; MIPGapAbs = abgap)
    solver = GurobiSolver(env)
    end

    #Bounds
    lb = [zeros(N,1); bnd[:,1]; zeros(d,1)]
    ub = [ones(N,1); bnd[:,2]; ones(d,1)]

    #Variable Type
    vtype = [repmat([:Bin],N); repmat([:Cont],k+d); repmat([:Bin],d)]

    if T > 0
        MathProgBase.setparameters!(solver, TimeLimit=T)
    end

    if abgap > 0
        setparams!(env;MIPGapAbs=abgap)
        solver = GurobiSolver(env)
    end

    ztemp1=zeros(N,d)
    ztemp2=zeros(2*d+1,N)
    htemp=[eye(d); -eye(d); zeros(1,d)]
    etemp=[-Diagonal(bnd[k+1:k+d,2]); Diagonal(bnd[k+1:k+d,1]); ones(1,d)]
    mtemp1=[ztemp2 zeros(2*d+1,k) htemp etemp]
    mtemp2=[zeros(2*d,1);[q]]

    if mio == 1 # Method 1 formulation
        obj = vcat((2*y-1), zeros(k+2*d,1))
        objcon = sum(1-y)
        miobnd_bar = miobnd+tol
        mtemp3= [-Diagonal(vec(miobnd_bar)) x_foc[:,2:k+1] x_aux ztemp1]
        A = [[Diagonal(vec(miobnd)) -x_foc[:,2:k+1] -x_aux ztemp1]; mtemp3; mtemp1]
        rhs = [miobnd*(1-tol)+beta0*x_foc[:,1];-tol*miobnd_bar-beta0*x_foc[:,1]; mtemp2]
    else # Method 2 formulation
        obj = [ones(1,N) zeros(1,k+2*d)]
        objcon = 0
        temp2=(1-2*y)
        A = [[Diagonal(vec(miobnd)) repmat(temp2,1,k).*x_foc[:,2:k+1] repmat(temp2,1,d).*x_aux ztemp1]; mtemp1]
        rhs = [miobnd*(1-tol)-(beta0*temp2.*x_foc[:,1]);mtemp2]
    end

    try
        m = buildmio(-vec(obj),A,'<',vec(rhs),vtype,vec(lb),vec(ub),solver)
        MathProgBase.optimize!(m)
        bhat = MathProgBase.getsolution(m)[N+1:N+k+d]
        score = -MathProgBase.getobjval(m) + objcon
        gap = (-MathProgBase.getobjbound(m)-score+objcon)
        rtime= MathProgBase.getsolvetime(m)
        ncount = MathProgBase.getnodecount(m)
        println("Optimization returned status: ", MathProgBase.status(m))

    catch
        println("Error reported")
    end
    return bhat[:,:], score, gap, rtime, ncount
end

export MaxScore


text:
module Tiger

import Base: print, isempty, convert, Iterators.flatten, Iterators.enumerate
import MAT.matread
using SparseArrays
using JuMP, Gurobi
using TOML

export parse_boolean,
       parse_gpr_string,
       read_cobra,
       build_base_model,
       read_media_file,
       set_media_bounds!,
       add_gprs_cnf!,
       single_deletions,
       extend_cobra_cnf,
       get_exchange_rxns,
       StandardLP,
       build_base_model

abstract type Boolean end
abstract type Junction <: Boolean end

struct Empty <: Boolean end

struct Atom <: Boolean
    name::String
end

struct And <: Junction
    left::Boolean
    right::Boolean
end

struct Or <: Junction
    left::Boolean
    right::Boolean
end

isempty(x::Empty) = true
isempty(x::Boolean) = false

convert(::Type{String}, a::Atom) = a.name

function print(x::Boolean)
    indent_step = "   "
    aux(a::Atom, indent) = println(indent * "\$" * a.name)
    aux(b::Or, indent) = begin
        println(indent * "Or")
        aux(b.left, indent * indent_step)
        aux(b.right, indent * indent_step)
    end
    aux(b::And, indent) = begin
        println(indent * "And")
        aux(b.left, indent * indent_step)
        aux(b.right, indent * indent_step)
    end
    aux(x, "")
end

function substitute_op(e)
    if isa(e, Expr) && e.head == :call
        if e.args[1] == :&
            e.args[1] = :And
        elseif e.args[1] == :|
            e.args[1] = :Or
        end

        for i = 2:length(e.args)
            e.args[i] = substitute_op(e.args[i])
        end
    elseif isa(e, Symbol) || isa(e, String)
        e = Atom(String(e))
    end

    return e
end

function parse_boolean(s::AbstractString)
    if isempty(s)
        return Empty()
    end
    ex = Meta.parse(s)
    return eval(substitute_op(ex))
end

insert(e::Union{Atom,And}, i::Or) = Or(i, e)
insert(e::Or, i::Or) = begin
    if isa(e.left, Or)
        insert(e.left, i)
    else
        Or(Or(i, e.left), e.right)
    end
end

dnf(e::Atom) = e
dnf(e::Or) = begin
    l = dnf(e.left)
    r = dnf(e.right)
    if isa(l, Atom) || isa(l, And)
        Or(r, l)
    elseif isa(r, Atom) || isa(r, And)
        Or(l, r)
    else
        insert(l, r)
    end
end
dnf(e::And) = begin
    l = dnf(e.left)
    r = dnf(e.right)
    if !isa(l, Or) && !isa(r, Or)
        e
    elseif isa(l, Or) && !isa(r, Or)
        Or(
            dnf(And(l.left, r)),
            And(l.right, r)
        )
    elseif isa(r, Or) && !isa(l, Or)
        Or(
            dnf(And(r.left, l)),
            And(r.right, l)
        )
    else
        # Or -- Or
        sub = Or(
            Or(
                And(l.left, r.right),
                And(l.left, r.left)
            ),
            And(l.right, r.left)
        )
        Or(
            dnf(sub),
            And(l.right, r.right)
        )
    end
end

swap_andor(b::Boolean) = b
swap_andor(b::And) = Or(swap_andor(b.left), swap_andor(b.right))
swap_andor(b::Or) = And(swap_andor(b.left), swap_andor(b.right))

cnf(b::Boolean) = b |> swap_andor |> dnf |> swap_andor

atoms(e::Empty) = []
atoms(a::Atom) = [a]
atoms(b::Boolean) = unique([atoms(b.left); atoms(b.right)])

name(a::Atom) = a.name
names(a::Atom) = [a.name]
names(b::Boolean) = name.(atoms(b))

right_names(b::Boolean) = [names(b)]
right_names(j::Junction) = push!(right_names(j.left), names(j.right))

dnf_groups(b::Boolean) = begin
    groups = []
    current = b
    while isa(current, Or)
        push!(groups, names(current.right))
        current = current.left
    end
    if !isempty(current)
        push!(groups, names(current))
    end
    return groups
end

cnf_groups(b::Boolean) = b |> swap_andor |> dnf_groups

# =================== Cobra Stuff ===================

function parse_gpr_string(s; and_str="and", or_str="or")
    if and_str != "&"
        s = replace(s, and_str => "&")
    end
    if or_str != "|"
        s = replace(s, or_str => "|")
    end
    parse_boolean(s)
end

struct CobraModel
    description::String

    c::Vector{Real}
    S::SparseMatrixCSC{Float64,Int}
    b::Vector{Real}
    lb::Vector{Real}
    ub::Vector{Real}
    csense::Vector{Char}

    rev::Vector{Bool}

    mets::Vector{String}
    metNames::Vector{String}
    rxns::Vector{String}
    rxnNames::Vector{String}

    vars::Vector{String}

    grRules::Vector{String}
    gprs::Vector{Boolean}
    genes::Vector{String}
end

"""
    read_cobra(matfile::AbstractString, name::AbstractString)

Read a Cobra model from a Matlab .mat file.

`name` is the variable name in the MAT file containing the model structure.
"""
function read_cobra(matfile::AbstractString, name::AbstractString)
    vars = matread(matfile)[name]
    gprs = parse_gpr_string.(vars["grRules"][:])
    genes = unique(flatten(atoms.(gprs)))
    CobraModel(
        vars["description"],

        vars["c"][:],
        sparse(vars["S"]),
        vars["b"][:],
        vars["lb"][:],
        vars["ub"][:],
        repeat(['='], length(vars["b"])),

        vars["rev"][:],

        vars["mets"][:],
        vars["metNames"][:],
        vars["rxns"][:],
        vars["rxnNames"][:],
        vars["rxns"][:],

        vars["grRules"][:],
        gprs,
        genes
    )
end

"""
    expand_cobra(cobra, ncons=0, nvars=0)

Add `ncons` constraints and `nvars` variables to a Cobra model.

Only the fields c, S, b, lb, ub, and csense are extended. If a variable is added,
all entries are zero. If a constraint is added, it take the form `0=0`.
"""
function expand_cobra(cobra; ncons=0, nvars=0, vars=nothing)
    m, n = size(cobra.S)

    c = zeros(n+nvars)
    c[1:n] = cobra.c

    S = spzeros(m+ncons, n+nvars)
    S[1:m,1:n] = cobra.S

    b = zeros(m+ncons)
    b[1:m] = cobra.b

    lb = zeros(n+nvars)
    lb[1:n] = cobra.lb

    ub = zeros(n+nvars)
    ub[1:n] = cobra.ub

    csense = repeat(['='], m+ncons)
    csense[1:m] = cobra.csense

    vars_ = repeat([""], n+nvars)
    vars_[1:n] = cobra.vars
    if !isnothing(vars)
        vars_[n+1:end] = vars
    end

    CobraModel(
        cobra.description,
        c,
        S,
        b,
        lb,
        ub,
        csense,
        cobra.rev,
        cobra.mets,
        cobra.metNames,
        cobra.rxns,
        cobra.rxnNames,
        vars_,
        cobra.grRules,
        cobra.gprs,
        cobra.genes
    )
end

"""
    build_base_model(cobra::CobraModel, optimizer=Gurobi.Optimizer, model=nothing)

Convert a CobraModel object into a JuMP model.

Adds variables for each reaction and defines the problem
    max c'*v
    s.t.
        S*v <|=|> b
        lb <= v <= ub

If model=nothing, a new model is created; otherwise the constraints are added to
the existing model.
"""
function build_base_model(cobra::CobraModel; optimizer=Gurobi.Optimizer, model=nothing)
    if isnothing(model)
        model = Model(optimizer)
    end

    nv = length(cobra.lb)

    #v = @variable(model, cobra.lb[i] <= [i = 1:nv] <= cobra.ub[i])
    v = @variable(model, [1:nv])
    set_name.(v, cobra.vars)
    set_lower_bound.(v, cobra.lb)
    set_upper_bound.(v, cobra.ub)

    sense = cobra.csense .== '<'
    if any(sense)
        @constraint(model, cobra.S[sense,:]*v .<= cobra.b[sense])
    end

    sense = cobra.csense .== '='
    if any(sense)
        @constraint(model, cobra.S[sense,:]*v .== cobra.b[sense])
    end

    sense = cobra.csense .== '>'
    if any(sense)
        @constraint(model, cobra.S[sense,:]*v .>= cobra.b[sense])
    end

    @objective(model, Max, cobra.c' * v)

    return model
end

"""
    add_gprs_cnf!(model::Model, cobra::CobraModel)

Add GPRs to a JuMP model using a CNF encoding.

Adds a continuous, nonnegative variable for each gene and constraints
to enforce the GPR rules. Variables are named based on the strings in
`cobra.genes`.

Setting a gene variable to zero is analogous to a gene knockout; otherwise
the value of each gene corresponds to the magnitude of the flux through the
the associated reactions. Note that the GPR logic is encoded. If two enzymes
can catalyze a reaction, the sum of the corresponding gene variables will
equal the reaction flux.
"""
function add_gprs_cnf!(model::Model, cobra::CobraModel)
    ng = length(cobra.genes)
    for gene in cobra.genes
        @variable(model, base_name=gene, lower_bound=0.0)
    end

    for (i, rxn) in enumerate(cobra.rxns)
        flux = variable_by_name(model, rxn)
        rule = cobra.gprs[i]
        if isempty(rule) continue end
        groups = cnf_groups(rule)
        for group in groups
            activities = variable_by_name.(model, group)
            @constraint(model, flux <= sum(activities))
            @constraint(model, flux >= -sum(activities))
        end
    end
end

function extend_cobra_cnf(cobra::CobraModel; ub=Inf)
    ng = length(cobra.genes)
    m, n = size(cobra.S)
    groups = cnf_groups.(cobra.gprs)
    ncons = 2(sum(length.(groups)))

    cobra = expand_cobra(cobra; ncons, nvars=ng, vars=cobra.genes)
    cobra.ub[n+1:end] .= ub

    gene_index(g) = findfirst(x -> x==g, cobra.vars)
    current = m
    for i = 1:n
        if isempty(groups[i]) continue end
        for group in groups[i]
            gidxs = gene_index.(group)
            current += 1
            cobra.S[current,i] = 1
            cobra.S[current,gidxs] .= -1
            cobra.csense[current] = '<'
            cobra.b[current] = 0

            current += 1
            cobra.S[current,i] = 1
            cobra.S[current,gidxs] .= 1
            cobra.csense[current] = '>'
            cobra.b[current] = 0
        end
    end

    return cobra
end

function read_media_file(mediafile::String)
    return TOML.parsefile(mediafile)
end

function set_media_bounds!(model::Model, media::Dict{String, Any}; set_defaults=true)
    if set_defaults
        varnames = all_variables(model) .|> JuMP.name
        for nm in varnames
            if occursin(media["exchange_pattern"], nm)
                var = variable_by_name(model, nm)
                set_lower_bound(var, media["default_exchange_bounds"][1])
                set_upper_bound(var, media["default_exchange_bounds"][2])
            end
        end
    end

    bounds = media["reactions"]
    for nm in keys(bounds)
        var = variable_by_name(model, nm)
        set_lower_bound(var, bounds[nm][1])
        set_upper_bound(var, bounds[nm][2])
    end
end

function set_media_bounds!(model::CobraModel, media::Dict{String, Any}; set_defaults=true)
    exchanges = get_exchange_rxns(model)
    for (i,rxn) in enumerate(model.rxns)
        if !(rxn in exchanges || haskey(media["reactions"], rxn))
            continue
        end

        if haskey(media["reactions"], rxn)
            model.lb[i] = media["reactions"][rxn][1]
            model.ub[i] = media["reactions"][rxn][2]
        elseif set_defaults
            model.lb[i] = media["default_exchange_bounds"][1]
            model.ub[i] = media["default_exchange_bounds"][2]
        end
    end
end

"""
    set_media_bounds!(model, filename::AbstractString; set_defaults=true)

Set exchange bounds using values in a TOML file.

See `./test/SampleMedia.toml` for a description of the media file format.
The media file can also be read using `read_media_file`. The resulting Dict
can be passed instead of the filename.
"""
function set_media_bounds!(model::Union{Model,CobraModel}, filename::AbstractString; set_defaults=true)
    set_media_bounds!(model, read_media_file(filename); set_defaults)
end

struct Bounds
    lb::Union{Nothing,Float64}
    ub::Union{Nothing,Float64}
end

"""
    save_bounds(var::VariableRef)

Create a Bounds object with `var`'s current bounds. The a bound is `nothing` if no bound is set.
"""
function save_bounds(var::VariableRef)
    Bounds(
        has_lower_bound(var) ? lower_bound(var) : nothing,
        has_upper_bound(var) ? upper_bound(var) : nothing
    )
end

"""
    restore_bounds!(var::VariableRef, bounds::Bounds)

Set the bounds of a JuMP variable. If either bound is `nothing`, the current bound is deleted.
"""
function restore_bounds!(var::VariableRef, bounds::Bounds)
    if isnothing(bounds.lb) && has_lower_bound(var)
        delete_lower_bound(var)
    else
        set_lower_bound(var, bounds.lb)
    end

    if isnothing(bounds.ub) && has_upper_bound(var)
        delete_upper_bound(var)
    else
        set_upper_bound(var, bounds.ub)
    end
end

"""
    single_deletions(model::Model, vars=all_variables(model))

Solve the model when each variable is constrainted to zero.

Returns vectors with the objective values for each deletion and the solver status.
"""
function single_deletions(model::Model, vars=all_variables(model))
    n = length(vars)
    objval = zeros(n)
    status = Array{Any}(undef, n)

    for (i, var) in enumerate(vars)
        prev_bounds = save_bounds(var)
        set_lower_bound(var, 0.0)
        set_upper_bound(var, 0.0)
        optimize!(model)
        objval[i] = objective_value(model)
        status[i] = termination_status(model)
        restore_bounds!(var, prev_bounds)
    end

    return objval, status
end

struct StandardLP
    # min sense * c'x
    # s.t.  Ax <= b
    #       Aeq x == beq
    #       lb <= x <= ub
    sense::Integer

    c::Vector{Real}

    A::SparseMatrixCSC{Float64,Int}
    b::Vector{Real}

    Aeq::SparseMatrixCSC{Float64,Int}
    beq::Vector{Real}

    lb::Vector{Real}
    ub::Vector{Real}

    vars::Vector{String}
end

function build_base_model(lp::StandardLP; optimizer=Gurobi.Optimizer, model=nothing)
    if isnothing(model)
        model = Model(optimizer)
    end

    n = length(lp.vars)
    x = @variable(model, [1:n])
    set_lower_bound.(x, lp.lb)
    set_upper_bound.(x, lp.ub)

    @constraint(model, lp.A*x .<= lp.b)
    @constraint(model, lp.Aeq*x .== lp.beq)

    if lp.sense == 1
        @objective(model, Min, lp.c'*x)
    else
        @objective(model, Max, lp.c'*x)
    end

    return model
end

function StandardLP(cobra::Tiger.CobraModel)
    cobra = deepcopy(cobra)
    sense = -1  # Cobra models maximize by default
    c = cobra.c
    lb = cobra.lb
    ub = cobra.ub
    n = size(cobra.S,2)

    ieq = cobra.csense .== '='
    Aeq = cobra.S[ieq,:]
    beq = cobra.b[ieq]

    leq = cobra.csense .== '<'
    Ale = cobra.S[leq,:]
    ble = cobra.b[leq]

    geq = cobra.csense .== '>'
    Age = cobra.S[geq,:]
    bge = cobra.b[geq]

    A = [Ale; -Age]
    b = [ble; -bge]

    StandardLP(
        sense,
        c,
        A,
        b,
        Aeq,
        beq,
        lb,
        ub,
        cobra.vars
    )
end

function get_exchange_rxns(cobra::Tiger.CobraModel; objective=true)
    # Extending a CobraModel, e.g. by adding genes via CNF, adds columns to the S matrix.
    # We only want to search in the original variables at indices 1:length(cobra.rxns).
    ex_idx = (sum(abs.(cobra.S[:,1:length(cobra.rxns)]),dims=1)  .== 1)[:]
    if !objective
        ex_idx .&= abs.(cobra.c[1:length(cobra.rxns)]) .== 0.0
    end
    cobra.rxns[ex_idx]
end

end  # module


text:
module CreateTableUsers

import SearchLight.Migrations: create_table, column, primary_key, add_index, drop_table

function up()
  create_table(:users) do
    [
      primary_key()
      column(:username, :string, limit = 100)
      column(:password, :string, limit = 100)
      column(:name, :string, limit = 100)
      column(:email, :string, limit = 100)
    ]
  end

  add_index(:users, :username)
end

function down()
  drop_table(:users)
end

end

max_stars_repo_name: genieframework/GenieAuthentication

max_issues_repo_name: genieframework/GenieAuthentication

max_forks_repo_name: genieframework/GenieAuthentication


text:
using FractionalDiffEq, Plots
plotlyjs()
(X, Y) = meshgrid(collect(-6:0.2:6), collect(-6:0.2:6))

Z = X.+Y.*im

L = mittleff(0.8, 0.9, Z)
L1=real.(L)
id = findall(x -> x>3, L1)
L1[id].=3
ix = findall(x-> x<-3, L1)
L1[ix].=-3
plot(X, Y, L1, st=:surface)


text:
# From Makie examples

@recipe(CorrPlot) do scene
    Theme(
        link = :x,  # need custom linking for y
        legend = false,
        margin = 1,
        fillcolor = theme(scene, :color),
        linecolor = theme(scene, :color),
        # indices = reshape(1:n^2, n, n)',
        title = "",
    )
end

AbstractPlotting.convert_arguments(::Type{<: CorrPlot}, x) = (x,)

function AbstractPlotting.plot!(scene::Scene, ::Type{CorrPlot}, attributes::Attributes, mat)
    n = size(mat, 2)
    layout = GridLayout(scene, n, n)
    layout[1:n, 1:n] = axs = [LAxis(scene) for i in 1:n, j in 1:n]

    for i in 1:n
        for j in 1:n
            ax = axs[i, j]
            if i > j
                scatter!(ax, view(mat, :, j), view(mat, :, i))
                plot!(ax, linear, view(mat, :, j), view(mat, :, i))
            elseif i == j
                plot!(ax, histogram, view(mat, :, j))
            else
                plot!(ax, histogram, view(mat, :, j), view(mat, :, i))
            end
        end
    end

    tight_ticklabel_spacing!.(axs)

    return scene
end


text:
function materialize!(curs::Cursor, restype::DataType)
    Co = typeof(curs).parameters[2]
    Fi = typeof(curs).parameters[3]
    n = ncols(curs)
    columns, nrows = _materialize(curs, restype)
    return columns, Fi, nrows
end


text:
module Delegates

export Delegator, Delegate

struct Delegator{T}
    delegate::T
end

struct Delegate end

operation(x::Delegator) = thing(x.delegate)
thing(::Any) = "default implementation"
thing(::Delegate) = "delegate implementation"

end  # module Delegates

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#
# Calculate molar weight from molecular formula
#

function molarweight(mform::String)
    reg = r"(?<atom>[a-zA-Z]+)(?<index>[0-9]+)"
    cs = 1
    m = match(reg, mform, cs)
    mw = 0.
    while m ≠ nothing
        aw = getaw(m[:atom])
        idx = parse(Int64, m[:index])
        mw += idx*aw
        cs += length(m.match)
        m = match(reg, mform, cs)
    end
    println("$mw g.mol⁻¹")
    return mw
end

function molarweight()
    mform = askmform()
    molarweight(mform)
end


text:
using Colors, CSV, DataFrames, PGFPlotsX
se_dmon = CSV.File("data/se_dmon1851_2019.tsv", delim="\t") |> DataFrame

function plotmonths(sex, months, syr, eyr)
    sax = @pgf Axis({"yticklabel style"={"/pgf/number format/use comma"},
              "xticklabel style"={"/pgf/number format/set thousands separator={}"},
              "yticklabel style"={"/pgf/number format/set thousands separator={}"},
              "legend style"={"font=\\tiny"},
              legend_pos="outer north east", ylabel="Döda",
              title = "Antal dödsfall per månad $(sex) Sverige",
              xmajorgrids, ymajorgrids})
    plotcolors = distinguishable_colors(length(months)+1, [RGB(1,1,1)])[2:end]
    for (i, month) in enumerate(months)
        plotmonth(sex, month, syr, eyr, sax, plotcolors[i])
    end
    sax
end

function plotmonth(sex, month, syr, eyr, p, plotcolor)
    se_dmon_sub = se_dmon[((se_dmon[:kön].==sex) .& (se_dmon[:månad].==month)
        .& (se_dmon[:år].>=syr) .& (se_dmon[:år].<=eyr)), :]
    @pgf push!(p, PlotInc({no_markers, color = plotcolor},
        Table([se_dmon_sub[:år], se_dmon_sub[:Döda]])),
        LegendEntry(month))
end

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Conversions related to the orbit elements.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Schwarz, R (2014). Memorandum No. 2: Cartesian State Vectors to
#       Keplerian Orbit Elements. Available at www.rene-schwarz.com
#
#       https://downloads.rene-schwarz.com/dc/category/18
#       (Accessed on 2017-08-09).
#
#   [2] Vallado, D. A., McClain, W. D (2013). Fundamentals of Astrodynamics
#       and Applications. Microcosm Press.
#
#   [3] Kuga, H. K., Carrara, V., Rao, K. R (2005). Introdução à Mecânica
#       Orbital. 2ª ed. Instituto Nacional de Pesquisas Espaciais.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export change_oe_frame

"""
    change_oe_frame(a::Number, e::Number, i::Number, Ω::Number, ω::Number, f::Number, conv_args...)
    change_oe_frame(oe::Orbit, conv_args...)

Change the reference frame of orbit elements. The orbit elements can be
specified by `a`, `e`, `i`, `Ω`, `ω`, and `f`, or the structure `oe` (see
[`Orbit`](@ref)). In the latter, the return value type will match the type of
`oe`.

The conversion arguments `conv_args` are **the same** arguments that one should
pass to the function [`r_eci_to_eci`](@ref) to convert between the desired
frames. For more information, see the documentation of the function
[`r_eci_to_eci`](@ref).

# Args

- `a`: Semi-major axis [m].
- `e`: Excentricity.
- `i`: Inclination [rad].
- `Ω`: Right-ascension of the ascending node [rad].
- `ω`: Argument of perigee [rad].
- `f`: True anomaly [rad].
- `conv_args...`: Conversion arguments, which are the same arguments that one
    would pass to the function `r_eci_to_eci` to convert between the desired
    frames.

- `oe`: An instance of the structure [`Orbit`](@ref) with the orbit elements
    that will be converted [SI units].

# Returns

Using the first signature, this function returns an instance of
[`KeplerianElements`](@ref). If the second signature is used, then the function
return an element with the same type of the input.

# Examples

```julia-repl
julia> eop = get_iers_eop(:IAU1980);

julia> teme_epoch = date_to_jd(2016,6,1,11,0,0);

julia> tod_epoch = date_to_jd(2016,1,1,0,0,0);

julia> k_teme = KeplerianElements(teme_epoch,
                                  7130.982e3,
                                  0.001111,
                                  98.405*pi/180,
                                  227.336*pi/180,
                                  90*pi/180,
                                  320*pi/180)
KeplerianElements{Float64}:
           Epoch :    2.45754e6 (2016-06-01T11:00:00)
 Semi-major axis : 7130.98     km
    Eccentricity :    0.001111
     Inclination :   98.405    °
            RAAN :  227.336    °
 Arg. of Perigee :   90.0      °
    True Anomaly :  320.0      °

julia> k_j2000 = change_oe_frame(k_teme, TEME(), J2000(), teme_epoch, eop)
KeplerianElements{Float64}:
           Epoch :    2.45754e6 (2016-06-01T11:00:00)
 Semi-major axis : 7130.98     km
    Eccentricity :    0.001111
     Inclination :   98.3365   °
            RAAN :  227.134    °
 Arg. of Perigee :   90.0604   °
    True Anomaly :  320.0      °

julia> k_tod = change_oe_frame(k_teme, TEME(), teme_epoch, TOD(), tod_epoch, eop)
KeplerianElements{Float64}:
           Epoch :    2.45754e6 (2016-06-01T11:00:00)
 Semi-major axis : 7130.98     km
    Eccentricity :    0.001111
     Inclination :   98.4037   °
            RAAN :  227.331    °
 Arg. of Perigee :   90.0014   °
    True Anomaly :  320.0      °
```
"""
function change_oe_frame(
        a::Number,
        e::Number,
        i::Number,
        Ω::Number,
        ω::Number,
        f::Number,
        conv_args...
)
    # The approach is to transform the orbit elements to Cartesian
    # representation, convert the frame, and then convert back to orbit
    # elements.

    k = KeplerianElements(0, a, e, i, Ω, ω, f)

    return change_oe_frame(k, conv_args...)
end

function change_oe_frame(oe::T, conv_args...) where T<:Orbit
    # First, we need to convert to state vector.
    sv_o = convert(OrbitStateVector, oe)
    r_o  = sv_o.r
    v_o  = sv_o.v
    a_o  = sv_o.a

    # NOTE: In my benchmarks, the operation with DCMs are faster than with
    # quaternions after the DCM representation was changed to SMatrix.
    D_ecif_ecio = r_eci_to_eci(DCM, conv_args...)
    r_f         = D_ecif_ecio * r_o
    v_f         = D_ecif_ecio * v_o
    a_f         = D_ecif_ecio * a_o

    # Create the new state vector that will be converted to an entity with the
    # type `T`.
    sv_f = orbsv(sv_o.t, r_f, v_f, a_f)

    return convert(T, sv_f)
end


text:
# https://oeis.org/A047662
# https://oeis.org/A008288
# https://oeis.org/A047665
# Fredman, "On the complexity of maintaining an array and its partial sums"
# tribonacci triangle (modified, because we have inhomogenous version)

using IndexedTables
using FileIO
# using LatexPrint
# using SurrealNumbers
out_dir = "Data/"

filename = "$(out_dir)multiplication_tables.csv"
t1 = load(filename)
t2 = table(t1)

# http://juliadb.org/latest/api/selection.html
# reindex(t2, (:n2))
# select(t2, :n2)
# filter(p -> eval(parse(p.n1)) < 3, t2)
# columns(t2, :n1)

n = 6
G = zeros(Int,n,n)
a = 0
i = 0
j = 0
G[1,:] = 1:n
G[:,1] = 1:n
for i=1:n
    for j=1:n
        a = filter(p -> eval(parse(p.n2)) == j, filter(p -> eval(parse(p.n1)) == i, t2))
        if !isempty(a)
            G[i,j] = columns(a, :generation)[1]
            G[j,i] = columns(a, :generation)[1]
        end
    end
end

m = 6
F = zeros(Int128,m,m)
F[1,:] = 1:m
F[:,1] = 1:m
for i=2:m
    for j=2:m
        F[i,j] = F[i-1,j] + F[i,j-1] + F[i-1,j-1] + 1
    end
end

N = zeros(Float64,m)
N2 = zeros(Float64,m)
λ = 3 + sqrt(2) # Mistake in Fredman
λ = 5.83 # approx
λ = 3 + 2*sqrt(2)
c = sqrt(λ/π) / 2^(9/4)
c2 = sqrt(8+6*sqrt(2))/(8*sqrt(π))
for k=1:m
    N[k] = c * λ^k / sqrt(k)
end

C = [ [A[k,k], N[k]] for k=1:m ]

# powers of 2
n = 0:6
theta = 1.597910218031873178338070118157
g_2_n = floor.(Int, theta.^(2.^n) )
# or f_n = f_n*(f_n + 1)
#    g(2^n) = f_n
# https://oeis.org/A007018
#


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Deprecation warnings.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

#                     Introduced in SatelliteToolbox v0.6
# ==============================================================================

@deprecate satellite_beta_angle(JD0, a, e, i, RAAN, numDays)                beta_angle(JD0, a, e, i, RAAN, numDays)
@deprecate satellite_check_station(r_e, rs_e, minElev)                      ground_station_visible(r_e, rs_e, minElev)
@deprecate satellite_check_station(r_e, lat_s, lon_s, h_s, minElev)         ground_station_visible(r_e, lat_s, lon_s, h_s, minElev)
@deprecate satellite_eclipse_time(JD0, a, e, i, w, RAAN, numDays)           eclipse_time_summary(JD0, a, e, i, RAAN, w, numDays)
@deprecate satellite_eclipse_time(JD0, a, e, i, w, RAAN, numDays, relative) eclipse_time_summary(JD0, a, e, i, RAAN, w, numDays, relative)
@deprecate igrf12(args...;kwargs...)                                        igrf(args...; kwargs...)

# Deprecation warnings related to Val{T} => Val(T)
# ==============================================================================

# Gravity models
# ------------------------------------------------------------------------------

@deprecate load_gravity_model(::Type{Val{T}}) where T load_gravity_model(Val(T))

# IGRF
# ------------------------------------------------------------------------------

@deprecate igrf12(date, r, λ, Ω, ::Type{Val{T}}; show_warns = true) where T igrf(date, r, λ, Ω, Val(T); show_warns = show_warns)

# Legendre
# ------------------------------------------------------------------------------

@deprecate dlegendre!(::Type{Val{T}}, dP, ϕ, P, ph_term = false) where T            dlegendre!(Val(T), P, ϕ, dP, ph_term)
@deprecate dlegendre(::Type{Val{T}}, ϕ, n_max, m_max = -1, ph_term = false) where T dlegendre(Val(T), ϕ, n_max, m_max, ph_term)
@deprecate legendre!(::Type{Val{T}}, P, ϕ, ph_term = false) where T                 legendre!(Val(T), P, ϕ, ph_term)
@deprecate legendre(::Type{Val{T}}, ϕ, n_max, m_max = -1, ph_term = false) where T  legendre(Val(T), ϕ, n_max, m_max, ph_term)

# Orbit propagation
# ------------------------------------------------------------------------------

# Individual parameters

@deprecate init_orbit_propagator(::Type{Val{:J2}}, epoch, a_0, e_0, i_0, Ω_0, ω_0, f_0, dn_o2 = 0, ddn_o6 = 0, j2_gc = j2_gc_egm08) where T #=
=#  init_orbit_propagator(Val(:J2), epoch, a_0, e_0, i_0, Ω_0, ω_0, f_0, dn_o2, ddn_o6, j2_gc)

@deprecate init_orbit_propagator(::Type{Val{:J4}}, epoch, a_0, e_0, i_0, Ω_0, ω_0, f_0, dn_o2 = 0, ddn_o6 = 0, j4_gc = j4_gc_egm08) where T #=
=#  init_orbit_propagator(Val(:J4), epoch, a_0, e_0, i_0, Ω_0, ω_0, f_0, dn_o2, ddn_o6, j4_gc)

@deprecate init_orbit_propagator(::Type{Val{:twobody}}, epoch, a_0, e_0, i_0, Ω_0, ω_0, f_0, μ = m0) where T #=
=#  init_orbit_propagator(Val(:twobody), epoch, a_0, e_0, i_0, Ω_0, ω_0, f_0, μ)

# Orbit structure

@deprecate init_orbit_propagator(::Type{Val{:J2}}, orb_0::Orbit, dn_o2::Number = 0, ddn_o6::Number = 0, j2_gc::J2_GravCte = j2_gc_egm08) #=
=#  init_orbit_propagator(Val(:J2), orb_0, dn_o2, ddn_o6, j2_gc)

@deprecate init_orbit_propagator(::Type{Val{:J4}}, orb_0::Orbit, dn_o2::Number = 0, ddn_o6::Number = 0, j4_gc::J4_GravCte = j4_gc_egm08) #=
=#  init_orbit_propagator(Val(:J4), orb_0, dn_o2, ddn_o6, j4_gc)

@deprecate init_orbit_propagator(::Type{Val{:twobody}}, orb_0::Orbit, μ::Number = m0) #=
=#  init_orbit_propagator(Val(:twobody), orb_0, μ)

# TLE

@deprecate init_orbit_propagator(::Type{Val{:J2}}, tle::TLE, j2_gc::J2_GravCte = j2_gc_egm08) #=
=#  init_orbit_propagator(Val(:J2), tle, j2_gc)

@deprecate init_orbit_propagator(::Type{Val{:J4}}, tle::TLE, j4_gc::J4_GravCte = j4_gc_egm08) #=
=#  init_orbit_propagator(Val(:J4), tle, j4_gc)

@deprecate init_orbit_propagator(::Type{Val{:twobody}}, tle::TLE, μ = m0) #=
=#  init_orbit_propagator(Val(:twobody), tle, μ)

@deprecate init_orbit_propagator(::Type{Val{:sgp4}}, tle::TLE, sgp4_gc::SGP4_GravCte = sgp4_gc_wgs84) #=
=#  init_orbit_propagator(Val(:sgp4), tle, sgp4_gc)

# Space indices
# ------------------------------------------------------------------------------

for sym in (:F10, :F10obs, :F10adj, :Kp, :Ap, :Kp_vect, :Ap_vect, :S10, :S81a,
            :M10, :M81a, :Y10, :Y81a, :DstΔTc)
    qsym = Meta.quot(sym)
    @eval @deprecate get_space_index(::Type{Val{$qsym}}, JD) get_space_index(Val($qsym), JD)
end

for sym in (:F10M, :F10Mobs, :F10Madj)
    qsym = Meta.quot(sym)
    @eval @deprecate get_space_index(::Type{Val{$qsym}}, JD; window = 81) get_space_index(Val($qsym), JD; window = window)
end

# Transformations
# ------------------------------------------------------------------------------

# ECEF to ECEF

@deprecate rECEFtoECEF(::Type{Val{To}}, ::Type{Val{Tf}}, JD, eop_data) where {To,Tf}    rECEFtoECEF(Val(To), Val(Tf), JD, eop_data)
@deprecate rECEFtoECEF(T, ::Type{Val{To}}, ::Type{Val{Tf}}, JD, eop_data) where {To,Tf} rECEFtoECEF(T, Val(To), Val(Tf), JD, eop_data)

# ECEF to ECI

@deprecate rECEFtoECI(::Type{Val{To}}, ::Type{Val{Tf}}, JD) where {To,Tf}              rECEFtoECI(Val(To), Val(Tf), JD)
@deprecate rECEFtoECI(::Type{Val{To}}, ::Type{Val{Tf}}, JD, eop_data) where {To,Tf}    rECEFtoECI(Val(To), Val(Tf), JD, eop_data)
@deprecate rECEFtoECI(T, ::Type{Val{To}}, ::Type{Val{Tf}}, JD) where {To,Tf}           rECEFtoECI(T, Val(To), Val(Tf), JD)
@deprecate rECEFtoECI(T, ::Type{Val{To}}, ::Type{Val{Tf}}, JD, eop_data) where {To,Tf} rECEFtoECI(T, Val(To), Val(Tf), JD, eop_data)

# ECI to ECEF

@deprecate rECItoECEF(::Type{Val{To}}, ::Type{Val{Tf}}, JD) where {To,Tf}              rECItoECEF(Val(To), Val(Tf), JD)
@deprecate rECItoECEF(::Type{Val{To}}, ::Type{Val{Tf}}, JD, eop_data) where {To,Tf}    rECItoECEF(Val(To), Val(Tf), JD, eop_data)
@deprecate rECItoECEF(T, ::Type{Val{To}}, ::Type{Val{Tf}}, JD) where {To,Tf}           rECItoECEF(T, Val(To), Val(Tf), JD)
@deprecate rECItoECEF(T, ::Type{Val{To}}, ::Type{Val{Tf}}, JD, eop_data) where {To,Tf} rECItoECEF(T, Val(To), Val(Tf), JD, eop_data)

# ECI to ECI

@deprecate rECItoECI(::Type{Val{To}}, ::Type{Val{Tf}}, JD) where {To,Tf}                    rECItoECI(Val(To), Val(Tf), JD)
@deprecate rECItoECI(::Type{Val{To}}, JDo, ::Type{Val{Tf}}, JDf) where {To,Tf}              rECItoECI(Val(To), JDo, Val(Tf), JDf)
@deprecate rECItoECI(::Type{Val{To}}, ::Type{Val{Tf}}, JD, eop_data) where {To,Tf}          rECItoECI(Val(To), Val(Tf), JD, eop_data)
@deprecate rECItoECI(::Type{Val{To}}, JDo, ::Type{Val{Tf}}, JDf, eop_data) where {To,Tf}    rECItoECI(Val(To), JDo, Val(Tf), JDf, eop_data)
@deprecate rECItoECI(T::T_ROT, ::Type{Val{To}}, ::Type{Val{Tf}}, JD) where {To,Tf}                 rECItoECI(T, Val(To), Val(Tf), JD)
@deprecate rECItoECI(T::T_ROT, ::Type{Val{To}}, JDo, ::Type{Val{Tf}}, JDf) where {To,Tf}           rECItoECI(T, Val(To), JDo, Val(Tf), JDf)
@deprecate rECItoECI(T::T_ROT, ::Type{Val{To}}, ::Type{Val{Tf}}, JD, eop_data) where {To,Tf}       rECItoECI(T, Val(To), Val(Tf), JD, eop_data)
@deprecate rECItoECI(T::T_ROT, ::Type{Val{To}}, JDo, ::Type{Val{Tf}}, JDf, eop_data) where {To,Tf} rECItoECI(T, Val(To), JDo, Val(Tf), JDf, eop_data)

# State vector, ECEF to ECI

@deprecate svECEFtoECI(sv, ::Type{Val{To}}, ::Type{Val{Tf}}, JD)           where {To,Tf} svECEFtoECI(sv, Val(To), Val(Tf), JD)
@deprecate svECEFtoECI(sv, ::Type{Val{To}}, ::Type{Val{Tf}}, JD, eop_data) where {To,Tf} svECEFtoECI(sv, Val(To), Val(Tf), JD, eop_data)

# State vector, ECI to ECEF

@deprecate svECItoECEF(sv, ::Type{Val{To}}, ::Type{Val{Tf}}, JD)           where {To,Tf} svECItoECEF(sv, Val(To), Val(Tf), JD)
@deprecate svECItoECEF(sv, ::Type{Val{To}}, ::Type{Val{Tf}}, JD, eop_data) where {To,Tf} svECItoECEF(sv, Val(To), Val(Tf), JD, eop_data)


text:
#    Problem :
#    *********
#    The BROYDN7D problem by Oren.
#
#    Source:
#    Ph.L. Toint,
#    "Some numerical results using a sparse matrix updating formula in
#    unconstrained optimization",
#    Mathematics of Computation, vol. 32(114), pp. 839-852, 1978.
#
#    See also Buckley#84
#
#    BROYDN7D.SIF classification OUR2-AN-V-0
#
#    Number of variables must be greater than 4
#
#    This problem has a reputation of issues.
#
# Daniel Henderson, 08/2021

f = (x) -> begin
	n = length(x)
	fx = abs(1-2x[2]+(3-2x[1])x[1])^(7/3) + abs(1-x[n-1]+(3-2x[n])x[n])^(7/3)
	for i in 1:Int(n/2)
		fx += abs(x[i]+x[i+Int(n/2)])^(7/3)
	end
	for i in 2:n-1
		fx += abs(1-x[i-1]-2*x[i+1]+(3-2*x[i])*x[i])^(7/3)
	end
	return fx
end

g! = (g, x) -> begin
	n      = length(x)
	first  = -2.0*x[1] + 1 + (3.0-2.0*x[1])x[1]
	last   = -x[n-1] + 1 + (3.0-2.0*x[n])x[n]
	g[1]   = 7/3*abs(first)^(4/3) * ( first > 0 ? 1 : -1) * (3.0-4x[1])
	g[2]   = -14/3*abs(first)^(4/3) * ( first > 0 ? 1 : -1)
	g[n-1] = -7/3*abs(last)^(4/3) * ( last > 0 ? 1 : -1)
	g[n]   = 7/3*abs(last)^(4/3) * ( last > 0 ? 1 : -1) * (3.0-4.0x[n])
	last   = x[1] + x[Int(n/2)+1]
	g[1]  += 7/3*abs(last)^(4/3) * ( last > 0 ? 1 : -1)
	g[Int(n/2)+1]  += 7/3*abs(last)^(4/3) * ( last > 0 ? 1 : -1)

	for i in 2:Int(n/2)
		first   = 1 - x[i-1] - 2x[i+1] + (3-2x[i])x[i]
		last    = x[i] + x[i + Int(n/2)]
		g[i-1] += -7.0/3*abs(first)^(4/3) * (first > 0 ? 1 : -1)
		g[i]   += 7.0/3*abs(first)^(4/3) * (first > 0 ? 1 : -1)*(3.0-4.0x[i]) + 7.0/3*abs(last)^(4/3) * (last > 0 ? 1 : -1)
		g[i+1] += -14/3*abs(first)^(4/3) * (first > 0 ? 1 : -1 )
		g[i+Int(n/2)] += 7.0/3*abs(last)^(4/3) * (last > 0 ? 1 : -1)
	end

	for i in Int(n/2)+1:n-1
		first   = 1 - x[i-1] - 2x[i+1] + (3-2x[i])x[i]
		g[i-1] += -7.0/3*abs(first)^(4/3) * (first > 0 ? 1 : -1)
		g[i]   += 7.0/3*abs(first)^(4/3) * (first > 0 ? 1 : -1)*(3.0-4.0x[i])
		g[i+1] += -14/3*abs(first)^(4/3) * (first > 0 ? 1 : -1 )
	end

	return g
end

fg! = (g, x) -> begin
	n      = length(x)
	first  = -2.0*x[1] + 1 + (3.0-2.0*x[1])x[1]
	last   = -x[n-1] + 1 + (3.0-2.0*x[n])x[n]
	fx     = abs(first)^(7/3) + abs(last)^(7/3)
	g[1]   = 7/3*abs(first)^(4/3) * ( first > 0 ? 1 : -1) * (3.0-4x[1])
	g[2]   = -14/3*abs(first)^(4/3) * ( first > 0 ? 1 : -1)
	g[n-1] = -7/3*abs(last)^(4/3) * ( last > 0 ? 1 : -1)
	g[n]   = 7/3*abs(last)^(4/3) * ( last > 0 ? 1 : -1) * (3.0-4.0x[n])
	last   = x[1] + x[Int(n/2)+1]
	fx    += abs(last)^(7/3)
	g[1]  += 7/3*abs(last)^(4/3) * ( last > 0 ? 1 : -1)
	g[Int(n/2)+1]  += 7/3*abs(last)^(4/3) * ( last > 0 ? 1 : -1)

	for i in 2:Int(n/2)
		first  = 1 - x[i-1] - 2x[i+1] + (3-2x[i])x[i]
		last   = x[i] + x[i + Int(n/2)]
		fx 	   += abs(first)^(7/3) + abs(last)^(7/3)
		g[i-1] += -7.0/3*abs(first)^(4/3) * (first > 0 ? 1 : -1)
		g[i]   += 7.0/3*abs(first)^(4/3) * (first > 0 ? 1 : -1)*(3.0-4.0x[i]) + 7.0/3*abs(last)^(4/3) * (last > 0 ? 1 : -1)
		g[i+1] += -14/3*abs(first)^(4/3) * (first > 0 ? 1 : -1 )
		g[i+Int(n/2)] += 7.0/3*abs(last)^(4/3) * (last > 0 ? 1 : -1)
	end

	for i in Int(n/2)+1:n-1
		first   = 1 - x[i-1] - 2x[i+1] + (3-2x[i])x[i]
		fx 	   += abs(first)^(7/3)
		g[i-1] += -7.0/3*abs(first)^(4/3) * (first > 0 ? 1 : -1)
		g[i]   += 7.0/3*abs(first)^(4/3) * (first > 0 ? 1 : -1)*(3.0-4.0x[i])
		g[i+1] += -14/3*abs(first)^(4/3) * (first > 0 ? 1 : -1 )
	end

	return fx, g
end

init = (n::Int=5000) -> begin
	mod(n, 2) > 0 && @warn "BROYDN7D: number of variables must be even"
	q = max(1, div(n, 2))
	n = 2q

    return n, ones(n)
end

@warn "TODO: Issue in CUTEst? See comment https://github.com/JuliaSmoothOptimizers/OptimizationProblems.jl/blob/main/src/broydn7d.jl#L50 "
# TestSet["BROYDN7D"] = UncProgram("BDQRTIC", f, g!, fg!, init)


text:
function solve_max_seat_id(input)
	max_seat_id = 0
	for pass in input
		seat_id = parse(Int16, replace(replace(pass, ('B', 'R') => '1'), ('F', 'L')=>'0'), base = 2)
		if seat_id > max_seat_id
			max_seat_id = seat_id
		end
	end
	max_seat_id
end

function solve_my_seat_id(input)
	my_seat_id = max_seat_id = 0
	min_seat_id = 1023
	for pass in input
		seat_id = parse(Int16, replace(replace(pass, ('B', 'R') => '1'), ('F', 'L')=>'0'), base = 2)
		if seat_id > max_seat_id
			max_seat_id = seat_id
		elseif seat_id < min_seat_id
			min_seat_id = seat_id
		end
		my_seat_id ⊻= seat_id
	end
	for i in 1:min_seat_id - 1
		my_seat_id ⊻= i
	end
	for i in max_seat_id + 1:1023
		my_seat_id ⊻= i
	end
	my_seat_id
end

@assert solve_max_seat_id(eachline("input/test/test_day5.txt")) == 820 ["test case 1 failed"]

println("First star: " * string(solve_max_seat_id(eachline("input/day5.txt"))))
println("Second star: " * string(solve_my_seat_id(eachline("input/day5.txt"))))


text:
using LinSpaces2
using Base.Test

# nextfloat-based start point in limited Float32 precision
x1,xN,N = 0.1f0, 1.0f0, 50
for i=1:1000000
    @test linspace2(x1,xN,N)[1] == x1
    @test linspace2(x1,xN,N)[N] == xN
    x1 = nextfloat(x1)
end

# Uniformly distributed endpoints
for i=1:1000000
    x1,xN,N = rand(), rand(), rand(2:10000)
    @test linspace2(x1,xN,N)[1] == x1
    @test linspace2(x1,xN,N)[N] == xN
end

# Log distributed endpoints
for i=1:1000000
    x1,xN,N = -log(rand()), -log(rand()), rand(2:10000)
    @test linspace2(x1,xN,N)[1] == x1
    @test linspace2(x1,xN,N)[N] == xN
end


text:
function plotbwh(b, w, h,  P, ttot)
    l = @layout [a b; c ]
    h1 = heatmap(b,  aspect_ratio=:equal, xlims = (0, P.nx), ylims = (0, P.nx), title=@sprintf("b - t=%3.2f", ttot) )
    h2 = heatmap(w,  aspect_ratio=:equal, xlims = (0, P.nx), ylims = (0, P.nx), title=@sprintf("w - t=%3.2f", ttot)  )
    h3 = heatmap(h,         aspect_ratio=:equal, xlims = (0, P.nx), ylims = (0, P.nx), title=@sprintf("h - t=%3.2f", ttot)  )
    display(plot(h1, h2, h3, layout = l))
end

max_stars_repo_name: jhardenberg/bwh

max_issues_repo_name: jhardenberg/bwh

max_forks_repo_name: jhardenberg/bwh


text:
function main()
	currentDir = pwd()
	mkpath(string(currentDir,"/preBuild/QBasis/Fermionic"))
	mkpath(string(currentDir,"/preBuild/Creation/Fermionic"))
	mkpath(string(currentDir,"/preBuild/FT/Fermionic"))

	return currentDir
end

main()

max_stars_repo_name: Gawatz/exactQmodule

max_issues_repo_name: Gawatz/exactQmodule

max_forks_repo_name: Gawatz/exactQmodule


text:
module TestCore

include("preamble.jl")

@testset begin
    @test ismutable([])
    @test ismutable(Dict(:a=>1))
    @test !ismutable(())
    @test !ismutable((a=1,))
    @test !ismutable(SVector(0))
end

end  # module


text:
using ITensorGaussianMPS
using LinearAlgebra
using ITensors

# Electrons

# Half filling
N = 50
Nf = N

@show N, Nf

# Maximum MPS link dimension
_maxlinkdim = 200

@show _maxlinkdim

# DMRG cutoff
_cutoff = 1e-8

# Hopping
t = 1.0

# Electron-electron on-site interaction
U = 1.0

@show t, U

# Make the free fermion Hamiltonian for the up spins
ampo_up = AutoMPO()
for n in 1:N-1
  ampo_up .+= -t, "Cdagup", n, "Cup", n+1
  ampo_up .+= -t, "Cdagup", n+1, "Cup", n
end

# Make the free fermion Hamiltonian for the down spins
ampo_dn = AutoMPO()
for n in 1:N-1
  ampo_dn .+= -t, "Cdagdn", n, "Cdn", n+1
  ampo_dn .+= -t, "Cdagdn", n+1, "Cdn", n
end

# Hopping Hamiltonian with 2*N spinless fermions,
# alternating up and down spins
h = hopping_hamiltonian(ampo_up, ampo_dn)

# Get the Slater determinant
Φ = slater_determinant_matrix(h, Nf)

# Create an MPS from the slater determinant.
# For now it only works without Sz conservation, this will be supported soon.
s = siteinds("Electron", N; conserve_qns = true, conserve_sz = false)
println("Making free fermion starting MPS")
@time ψ0 = slater_determinant_to_mps(s, Φ; eigval_cutoff = 1e-4,
                                           cutoff = _cutoff,
                                           maxdim = _maxlinkdim)
@show maxlinkdim(ψ0)

@show U
ampo = ampo_up + ampo_dn
for n in 1:N
  ampo .+= U, "Nupdn", n
end
H = MPO(ampo, s)

# Random starting state
ψr = randomMPS(s, n -> n ≤ Nf ? (isodd(n) ? "↑" : "↓") : "0")

println("Random starting state energy")
@show flux(ψr)
@show inner(ψr, H, ψr)
println()
println("Free fermion starting state energy")
@show flux(ψ0)
@show inner(ψ0, H, ψ0)

println("\nStart from product state")
sweeps = Sweeps(10)
maxdim!(sweeps,10,20,_maxlinkdim)
cutoff!(sweeps,_cutoff)
@time dmrg(H, ψr, sweeps)

println("\nStart from free fermion state")
sweeps = Sweeps(5)
maxdim!(sweeps,_maxlinkdim)
cutoff!(sweeps,_cutoff)
@time dmrg(H, ψ0, sweeps)

nothing


text:
function cuMPO(O::MPO)
    P = copy(O)
    for site in 1:length(O)
        P.data[site] = cuITensor(O.data[site])
    end
    return P
end
cuMPO() = MPO()

cuMPO(A::Vector{ITensor}) = cuMPO(MPO(A))
cuMPO(sites) = cuMPO(MPO(sites))
function plussers(::Type{T}, left_ind::Index, right_ind::Index, sum_ind::Index) where {T <: CuArray}
    #if dir(left_ind) == dir(right_ind) == Neither
        total_dim    = dim(left_ind) + dim(right_ind)
        total_dim    = max(total_dim, 1)
        left_data   = CuArrays.zeros(Float64, dim(left_ind)*dim(sum_ind))
        ldi = diagind(reshape(left_data, dim(left_ind), dim(sum_ind)), 0)
        left_data[ldi] = 1.0
        left_tensor = cuITensor(vec(left_data), left_ind, sum_ind)
        right_data   = CuArrays.zeros(Float64, dim(right_ind) * dim(sum_ind))
        rdi = diagind(reshape(right_data, dim(right_ind), dim(sum_ind)), dim(left_ind))
        right_data[rdi] = 1.0
        right_tensor = cuITensor(vec(right_data), right_ind, sum_ind)
        return left_tensor, right_tensor
    #else # tensors have QNs
    #    throw(ArgumentError("support for adding MPOs with defined quantum numbers not implemented yet."))
    #end
end

function randomCuMPO(sites, m::Int=1)
  M = cuMPO(sites)
  for i ∈ eachindex(sites)
    randn!(M[i])
    normalize!(M[i])
  end
  m > 1 && throw(ArgumentError("randomMPO: currently only m==1 supported"))
  return M
end

function Base.collect(M::T) where {T <: Union{MPS, MPO}}
    if typeof(tensor(ITensors.data(M)[1])) <: CuDenseTensor
        return T(collect.(ITensors.data(M)), M.llim, M.rlim)
    else
        return M
    end
end