text:
using Test
using Pluto.Configuration: CompilerOptions
using Pluto.WorkspaceManager: _merge_notebook_compiler_options
import Pluto: update_save_run!, WorkspaceManager, ClientSession, ServerSession, Notebook, Cell, project_relative_path

@testset "Workspace manager" begin
# basic functionality is already tested by the reactivity tests

    @testset "Multiple notebooks" begin

        fakeclientA = ClientSession(:fakeA, nothing)
        fakeclientB = ClientSession(:fakeB, nothing)
        🍭 = ServerSession()
        🍭.connected_clients[fakeclientA.id] = fakeclientA
        🍭.connected_clients[fakeclientB.id] = fakeclientB

        notebookA = Notebook([
            Cell("x = 3")
        ])
        fakeclientA.connected_notebook = notebookA

        notebookB = Notebook([
            Cell("x")
        ])
        fakeclientB.connected_notebook = notebookB

        @test notebookA.path != notebookB.path

        update_save_run!(🍭, notebookA, notebookA.cells[1])
        update_save_run!(🍭, notebookB, notebookB.cells[1])

        @test notebookB.cells[1].errored == true
    end

    @testset "notebook environment" begin
        session_options = CompilerOptions()
        notebook = Notebook([Cell("x")])
        notebook.compiler_options = CompilerOptions(;project="test")
        @test _merge_notebook_compiler_options(notebook, session_options).project ==
            joinpath(dirname(notebook.path), "test")

        notebook.compiler_options = CompilerOptions(;project=project_relative_path("test"))
        @test _merge_notebook_compiler_options(notebook, session_options).project ==
            project_relative_path("test")

        session_options = CompilerOptions(;project=project_relative_path("test"))
        notebook.compiler_options = CompilerOptions(;project=project_relative_path("Project.toml"))
        @test _merge_notebook_compiler_options(notebook, session_options).project ==
            project_relative_path("Project.toml")
    end
end


text:
include("../common.jl")

# Offset to add to the version number
version_offset = v"0.1.0"

# Minimum Julia version supported: this is important to decide which versions of
# the dependencies to use, in particular the JLL stdlibs.
min_julia_version = v"1.3"

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, configure(version_offset, min_julia_version)...; julia_compat="~1.0, ~1.1, ~1.2, ~1.3, ~1.4, ~1.5")

ext: jl

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
include("../utils.jl")

using Flux, Flux.Tracker

const loop_iters = 500

function loop(x = 2, n = loop_iters)
  r = one(x)/one(x)
  while n > 0
    r *= sin(x)
    n -= 1
  end
  return r
end

@result loop Tracker gradient(loop, 2)

logreg(W, b, x) = sum(σ.(W * x .+ b))

W = randn(Float32, 10,28^2)
b = randn(Float32, 10)
x = rand(Float32, 28^2)

@result logreg Tracker gradient(logreg, $W, $b, $x)


text:
# Kernels for unary array operations

fp = open("cuda1.cu","w")
#using Knet: unary_ops

include("gamma.jl")
print(fp,cuda1gammafamily())

function cuda1src(f, j=f, ex="$f(xi)"; seperate_impl=false, BLK=256, THR=256)
    sprint() do s
        for (T,F) in [("float","$(f)_32"),("double","$(f)_64")]
            seperate_impl && (ex = "$F(xi)")
            print(s,
"""
__global__ void _$F(int n, $T *x, $T *y) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    $T xi = x[i];
    y[i] = $ex;
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  $DLLEXPORT void $F(int n, $T *x, $T *y) {
    if (n>0) _$F<<<$BLK,$THR>>>(n,x,y);
  }
}
""")
        end
    end
end

seperate_impl_ops = ["gamma_impl", "digamma_impl", "trigamma_impl"]
for a in unary_ops
    if !isa(a,Tuple); a=(a,); end
    print(fp,cuda1src(a...; seperate_impl=(a[1] in seperate_impl_ops)))
end

# Kernels used by setindex! and getindex: fill, xfill, xcopy:

function cuda1fill(; BLK=256, THR=256)
    sprint() do s
        for (T,F) in [("float","32"),("double","64")]
            print(s,
"""
__global__ void _fill_$F(int n, $T x, $T *y) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    y[i] = x;
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  $DLLEXPORT void fill_$F(int n, $T x, $T *y) {
    if (n>0) _fill_$F<<<$BLK,$THR>>>(n,x,y);
  }
}
""")
        end
    end
end

print(fp,cuda1fill())

function cuda1xfill(; BLK=256, THR=256)
    sprint() do s
        for (T,F) in [("float","32"),("double","64")]
            print(s,
"""
__global__ void _xfill_$F(int nrows, int ncols, $T x, $T *y, int incy) {
  int row, col, yidx;
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = i % nrows;
    col = i / nrows;
    if (col >= ncols) break;
    yidx = row + col * incy;
    y[yidx] = x;
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  $DLLEXPORT void xfill_$F(int nrows, int ncols, $T x, $T *y, int incy) {
    if (nrows>0 && ncols>0) _xfill_$F<<<$BLK,$THR>>>(nrows, ncols, x, y, incy);
  }
}
""")
        end
    end
end

print(fp,cuda1xfill())

function cuda1xcopy(; BLK=256, THR=256)
"""
__global__ void _xcopy(int nrows, int ncols, const char *x, int incx, char *y, int incy) {
  int row, col, xidx, yidx;
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = i % nrows;
    col = i / nrows;
    if (col >= ncols) break;
    xidx = row + col * incx;
    yidx = row + col * incy;
    y[yidx] = x[xidx];
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  $DLLEXPORT void xcopy(int nrows, int ncols, const void *x, int incx, void *y, int incy) {
    if (nrows>0 && ncols>0) _xcopy<<<$BLK,$THR>>>(nrows,ncols,(char*)x,incx,(char*)y,incy);
  }
}
"""
end

print(fp,cuda1xcopy())

### Kernels for permutedims by Ekrem Emre Yurdakul 2017-02-27

function permutedims2Dsrc(f,i1,i2; BLK=256,THR=256)
    sprint() do s
        for (T,F) in [("float","$(f)32"),("double","$(f)64")]
            print(s,
"""
__global__ void _$(F)($T* x, int dimx1, int dimx2, $T* y, int dimy1) {
  for (int v = threadIdx.x + blockIdx.x * blockDim.x; v < dimx1*dimx2; v += blockDim.x * gridDim.x) {

    //From 1D to 2D indices
    int i = v % dimx1;
    int j = (v-i) / dimx1;

    //Calculate destination
    int destIndex = $i1 + $i2*dimy1;
    y[destIndex] = x[v];
	}
}
extern "C" {
  $DLLEXPORT void $(F)($T* x, int dimx1, int dimx2, $T* y, int dimy1) {
    _$(F)<<<$BLK,$THR>>>(x,dimx1,dimx2,y,dimy1);
  }
}
""")
        end
    end
end

function permutedims3Dsrc(f,i1,i2,i3; BLK=256,THR=256)
    sprint() do s
        for (T,F) in [("float","$(f)32"),("double","$(f)64")]
            print(s,
"""
__global__ void _$(F)($T* x, int dimx1, int dimx2, int dimx3, $T* y, int dimy1, int dimy2) {
  for (int v = threadIdx.x + blockIdx.x * blockDim.x; v < dimx1*dimx2*dimx3; v += blockDim.x * gridDim.x) {

    //From 1D to 3D indices
    int i = v % dimx1;
    int j = ((v-i) / dimx1) % dimx2;
    int k = (v-i-j*dimx1) / (dimx1*dimx2);

    //Calculate destination
    int destIndex = $i1 + $i2*dimy1 + $i3*dimy1*dimy2;
    y[destIndex] = x[v];
	}
}
extern "C" {
  $DLLEXPORT void $(F)($T* x, int dimx1, int dimx2, int dimx3, $T* y, int dimy1, int dimy2) {
    _$(F)<<<$BLK,$THR>>>(x,dimx1,dimx2,dimx3,y,dimy1,dimy2);
  }
}
""")
        end
    end
end

function permutedims4Dsrc(f,i1,i2,i3,i4; BLK=256,THR=256)
    sprint() do s
        for (T,F) in [("float","$(f)32"),("double","$(f)64")]
            print(s,
"""
__global__ void _$(F)($T* x, int dimx1, int dimx2, int dimx3, int dimx4, $T* y, int dimy1, int dimy2, int dimy3) {
  for (int v = threadIdx.x + blockIdx.x * blockDim.x; v < dimx1*dimx2*dimx3*dimx4; v += blockDim.x * gridDim.x) {

    //From 1D to 4D indices
    int i = v % dimx1;
    int j = ((v-i) / dimx1) % dimx2;
    int k = ((v-i-j*dimx1) / (dimx1*dimx2)) % dimx3;
    int l = (v-i-j*dimx1-k*dimx1*dimx2) / (dimx1*dimx2*dimx3);

    //Calculate destination
    int destIndex = $i1 + $i2*dimy1 + $i3*dimy1*dimy2 + $i4*dimy1*dimy2*dimy3;
    y[destIndex] = x[v];
	}
}
extern "C" {
  $DLLEXPORT void $(F)($T* x, int dimx1, int dimx2, int dimx3, int dimx4, $T* y, int dimy1, int dimy2, int dimy3) {
    _$(F)<<<$BLK,$THR>>>(x,dimx1,dimx2,dimx3,dimx4,y,dimy1,dimy2,dimy3);
  }
}
""")
        end
    end
end

function permutedims5Dsrc(f,i1,i2,i3,i4,i5; BLK=256,THR=256)
    sprint() do s
        for (T,F) in [("float","$(f)32"),("double","$(f)64")]
            print(s,
"""
__global__ void _$(F)($T* x, int dimx1, int dimx2, int dimx3, int dimx4, int dimx5, $T* y, int dimy1, int dimy2, int dimy3, int dimy4) {
  for (int v = threadIdx.x + blockIdx.x * blockDim.x; v < dimx1*dimx2*dimx3*dimx4*dimx5; v += blockDim.x * gridDim.x) {

    //From 1D to 5D indices
    int i = v % dimx1;
    int j = ((v-i) / dimx1) % dimx2;
    int k = ((v-i-j*dimx1) / (dimx1*dimx2)) % dimx3;
    int l = ((v-i-j*dimx1-k*dimx1*dimx2) / (dimx1*dimx2*dimx3)) % dimx4;
    int m = (v-i-j*dimx1-k*dimx1*dimx2-l*dimx1*dimx2*dimx3) / (dimx1*dimx2*dimx3*dimx4);

    //Calculate destination
    int destIndex = $i1 + $i2*dimy1 + $i3*dimy1*dimy2 + $i4*dimy1*dimy2*dimy3 + $i5*dimy1*dimy2*dimy3*dimy4;
    y[destIndex] = x[v];
	}
}
extern "C" {
  $DLLEXPORT void $(F)($T* x, int dimx1, int dimx2, int dimx3, int dimx4, int dimx5, $T* y, int dimy1, int dimy2, int dimy3, int dimy4) {
    _$(F)<<<$BLK,$THR>>>(x,dimx1,dimx2,dimx3,dimx4,dimx5,y,dimy1,dimy2,dimy3,dimy4);
  }
}
""")
        end
    end
end

# avoid dependency, copy permutations from julia4 base
# using Combinatorics

import Base: iterate, length
struct Perms; a; end
perms(a) = Perms(a)
length(p::Perms) = factorial(length(p.a))
#start(p::Perms) = [1:length(p.a);]
#done(p::Perms, s) = !isempty(s) && s[1] > length(p.a)
function iterate(p::Perms, s=[1:length(p.a);])
    if !isempty(s) && s[1] > length(p.a); return nothing; end
    perm = [p.a[si] for si in s]
    if length(p.a) == 0
        # special case to generate 1 result for len==0
        return (perm,[1])
    end
    s = copy(s)
    k = length(s)-1
    while k > 0 && s[k] > s[k+1];  k -= 1;  end
    if k == 0
        s[1] = length(s)+1   # done
    else
        l = length(s)
        while s[k] >= s[l];  l -= 1;  end
        s[k],s[l] = s[l],s[k]
        reverse!(s,k+1)
    end
    (perm,s)
end

function cuda1permutedims()
  cudaPerms = [permutedims2Dsrc,permutedims3Dsrc,permutedims4Dsrc,permutedims5Dsrc]
  for i=2:5
      dims = collect(perms(1:i))
      indnames = collect(perms(["i","j","k","l","m"][1:i]))
      for j=1:length(dims)
          fname = string("permutedims_",i,"D_",join(dims[j],"_"),"_")
          print(fp,cudaPerms[i-1](fname,indnames[j]...))
      end
  end
end

cuda1permutedims()

function cuda1icat(; BLK=256, THR=256)
    sprint() do s
        for (T,F) in [("float","32"),("double","64")]
            print(s,
"""
__global__ void _icat_$F(int nrows, int ncols, $T **x, $T *y) {
  int row, col, yidx;
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = i % nrows;
    col = i / nrows;
    if (col >= ncols) break;
    yidx = row + col * nrows;
    y[yidx] = x[col][row];
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  $DLLEXPORT void icat_$F(int nrows, int ncols, $T **x, $T *y) {
    $T **xx;
    if (nrows>0 && ncols>0) {
      size_t s = ncols * sizeof($T *);
      cudaMalloc(&xx, s);
      cudaMemcpy(xx, x, s, cudaMemcpyHostToDevice);
      _icat_$F<<<$BLK,$THR>>>(nrows, ncols, xx, y);
      cudaFree(xx);
    }
  }
}
""")
        end
    end
end

print(fp,cuda1icat())

# This is for missing double atomicAdd()
print(fp,"""
static __inline__ __device__ float atomicAdd2(float *address, float val) {
  return atomicAdd(address, val);
}
#if defined(__CUDA_ARCH__) && __CUDA_ARCH__ >= 600
static __inline__ __device__ double atomicAdd2(double *address, double val) {
  return atomicAdd(address, val);
}
#else
static __inline__ __device__ double atomicAdd2(double *address, double val) {
  unsigned long long int* address_as_ull = (unsigned long long int*)address;
  unsigned long long int old = *address_as_ull, assumed;
  if (val==0.0)
    return __longlong_as_double(old);
  do {
    assumed = old;
    old = atomicCAS(address_as_ull, assumed, __double_as_longlong(val +__longlong_as_double(assumed)));
  } while (assumed != old);
  return __longlong_as_double(old);
}
#endif
""")

function cuda1getcols(; BLK=256, THR=256)
    sprint() do s
        for (T,F) in [("float","32"),("double","64")]
            print(s,
"""
__global__ void _getcols_$F(int xrows, int xcols, int ncols, int *cols, $T *x, $T *y) {
  int row, col, xidx;
  int yidx = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = yidx % xrows;
    col = yidx / xrows;
    if (col >= ncols) break;
    xidx = row + (cols[col]-1) * xrows;
    y[yidx] = x[xidx];
    yidx += blockDim.x * gridDim.x;
  }
}
__global__ void _setcols_$F(int xrows, int xcols, int ncols, int *cols, $T *x, $T *y) {
  int row, col, xidx;
  int yidx = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = yidx % xrows;
    col = yidx / xrows;
    if (col >= ncols) break;
    xidx = row + (cols[col]-1) * xrows;
    x[xidx] = y[yidx];
    yidx += blockDim.x * gridDim.x;
  }
}
__global__ void _addcols_$F(int xrows, int xcols, int ncols, int *cols, $T *x, $T *y) {
  int row, col, xidx;
  int yidx = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = yidx % xrows;
    col = yidx / xrows;
    if (col >= ncols) break;
    xidx = row + (cols[col]-1) * xrows;
    atomicAdd2(&x[xidx], y[yidx]);
    yidx += blockDim.x * gridDim.x;
  }
}
__global__ void _setcol1_$F(int xrows, int xcols, int ncols, int *cols, $T *x, $T y) {
  int row, col, xidx;
  int yidx = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = yidx % xrows;
    col = yidx / xrows;
    if (col >= ncols) break;
    xidx = row + (cols[col]-1) * xrows;
    x[xidx] = y;
    yidx += blockDim.x * gridDim.x;
  }
}
__global__ void _getrows_$F(int xrows, int xcols, int nrows, int *rows, $T *x, $T *y) {
  int row, col, xidx;
  int yidx = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = yidx % nrows;
    col = yidx / nrows;
    if (col >= xcols) break;
    xidx = rows[row] - 1 + col * xrows;
    y[yidx] = x[xidx];
    yidx += blockDim.x * gridDim.x;
  }
}
__global__ void _setrows_$F(int xrows, int xcols, int nrows, int *rows, $T *x, $T *y) {
  int row, col, xidx;
  int yidx = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = yidx % nrows;
    col = yidx / nrows;
    if (col >= xcols) break;
    xidx = rows[row] - 1 + col * xrows;
    x[xidx] = y[yidx];
    yidx += blockDim.x * gridDim.x;
  }
}
__global__ void _addrows_$F(int xrows, int xcols, int nrows, int *rows, $T *x, $T *y) {
  int row, col, xidx;
  int yidx = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = yidx % nrows;
    col = yidx / nrows;
    if (col >= xcols) break;
    xidx = rows[row] - 1 + col * xrows;
    atomicAdd2(&x[xidx], y[yidx]);
    yidx += blockDim.x * gridDim.x;
  }
}
__global__ void _setrow1_$F(int xrows, int xcols, int nrows, int *rows, $T *x, $T y) {
  int row, col, xidx;
  int yidx = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = yidx % nrows;
    col = yidx / nrows;
    if (col >= xcols) break;
    xidx = rows[row] - 1 + col * xrows;
    x[xidx] = y;
    yidx += blockDim.x * gridDim.x;
  }
}
__global__ void _getents_$F(int n, int *ents, $T *x, $T *y) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    y[i] = x[ents[i]-1];
    i += blockDim.x * gridDim.x;
  }
}
__global__ void _setents_$F(int n, int *ents, $T *x, $T *y) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    x[ents[i]-1] = y[i];
    i += blockDim.x * gridDim.x;
  }
}
__global__ void _addents_$F(int n, int *ents, $T *x, $T *y) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    atomicAdd2(&x[ents[i]-1], y[i]);
    i += blockDim.x * gridDim.x;
  }
}
__global__ void _setent1_$F(int n, int *ents, $T *x, $T y) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    x[ents[i]-1] = y;
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
$DLLEXPORT void getcols_$F(int xrows, int xcols, int ncols, int *cols, $T *x, $T *y)
{ if (ncols>0 && xrows>0 && xcols>0) _getcols_$F<<<$BLK,$THR>>>(xrows,xcols,ncols,cols,x,y); }
$DLLEXPORT void setcols_$F(int xrows, int xcols, int ncols, int *cols, $T *x, $T *y)
{ if (ncols>0 && xrows>0 && xcols>0) _setcols_$F<<<$BLK,$THR>>>(xrows,xcols,ncols,cols,x,y); }
$DLLEXPORT void addcols_$F(int xrows, int xcols, int ncols, int *cols, $T *x, $T *y)
{ if (ncols>0 && xrows>0 && xcols>0) _addcols_$F<<<$BLK,$THR>>>(xrows,xcols,ncols,cols,x,y); }
$DLLEXPORT void setcol1_$F(int xrows, int xcols, int ncols, int *cols, $T *x, $T  y)
{ if (ncols>0 && xrows>0 && xcols>0) _setcol1_$F<<<$BLK,$THR>>>(xrows,xcols,ncols,cols,x,y); }
$DLLEXPORT void getrows_$F(int xrows, int xcols, int nrows, int *rows, $T *x, $T *y)
{ if (nrows>0 && xrows>0 && xcols>0) _getrows_$F<<<$BLK,$THR>>>(xrows,xcols,nrows,rows,x,y); }
$DLLEXPORT void setrows_$F(int xrows, int xcols, int nrows, int *rows, $T *x, $T *y)
{ if (nrows>0 && xrows>0 && xcols>0) _setrows_$F<<<$BLK,$THR>>>(xrows,xcols,nrows,rows,x,y); }
$DLLEXPORT void addrows_$F(int xrows, int xcols, int nrows, int *rows, $T *x, $T *y)
{ if (nrows>0 && xrows>0 && xcols>0) _addrows_$F<<<$BLK,$THR>>>(xrows,xcols,nrows,rows,x,y); }
$DLLEXPORT void setrow1_$F(int xrows, int xcols, int nrows, int *rows, $T *x, $T  y)
{ if (nrows>0 && xrows>0 && xcols>0) _setrow1_$F<<<$BLK,$THR>>>(xrows,xcols,nrows,rows,x,y); }
$DLLEXPORT void getents_$F(int n, int *ents, $T *x, $T *y)
{ if (n>0) _getents_$F<<<$BLK,$THR>>>(n,ents,x,y); }
$DLLEXPORT void setents_$F(int n, int *ents, $T *x, $T *y)
{ if (n>0) _setents_$F<<<$BLK,$THR>>>(n,ents,x,y); }
$DLLEXPORT void addents_$F(int n, int *ents, $T *x, $T *y)
{ if (n>0) _addents_$F<<<$BLK,$THR>>>(n,ents,x,y); }
$DLLEXPORT void setent1_$F(int n, int *ents, $T *x, $T  y)
{ if (n>0) _setent1_$F<<<$BLK,$THR>>>(n,ents,x,y); }
}
""")
        end
    end
end

print(fp,cuda1getcols())

# Dropout

function cuda1dropout(; BLK=256, THR=256)
    sprint() do s
        for (T,F) in [("float","32"),("double","64")]
            print(s,
"""
__global__ void _dropout_$F(int n, $T p, $T q, $T *x, $T *y) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    if (y[i] > p) {
      y[i] = x[i] * q;
    } else {
      y[i] = 0;
    }
    i += blockDim.x * gridDim.x;
  }
}
__global__ void _dropback_$F(int n, $T q, $T *y, $T *dy, $T *dx) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    if (y[i] == 0) {
        dx[i] = 0;
    } else {
        dx[i] = dy[i] * q;
    }
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  $DLLEXPORT void dropout_$F(int n, $T p, $T *x, $T *y) {
    if (n>0) _dropout_$F<<<$BLK,$THR>>>(n,p,1.0/(1.0-p),x,y);
  }
  $DLLEXPORT void dropback_$F(int n, $T p, $T *x, $T *y, $T *dy, $T *dx) {
    if (n>0) _dropback_$F<<<$BLK,$THR>>>(n,1.0/(1.0-p),y,dy,dx);
  }
}
""")
        end
    end
end

print(fp,cuda1dropout())

# This is still too slow compared to concat on cpu and copy to gpu
# Tested for 25 arrays of 200
function cuda1concat(; BLK=256, THR=256)
    sprint() do s
        for (T,F) in [("float","32"),("double","64")]
            print(s,
"""
__global__ void _concat_$F(int narrays, int *starts, int *lengths, $T **x, $T *y) {
  int array = blockIdx.x;
  int nelts = lengths[array];
  int offset = starts[array];
  for (int i = threadIdx.x; i < nelts; i += blockDim.x) {
    y[i+offset] = x[array][i];
  }
}
extern "C" {
  // julia is responsible for copying args to gpu
  $DLLEXPORT void concat_$F(int narrays, int *starts, int *lengths, $T **x, $T *y) {
    _concat_$F<<<narrays,$THR>>>(narrays, starts, lengths, x, y);
  }
}
""")
        end
    end
end

print(fp,cuda1concat())

close(fp)

# Here is the test script for cuda1concat:
# using Knet, BenchmarkTools
# using Knet: @knet8

# for S in (32,64); T = Symbol("Float$S"); F = "concat_$S"
# @eval function concat(A::KnetArray{$T}...)
#     nargs = length(A)
#     S = Array{Int32}(nargs)
#     L = Array{Int32}(nargs)
#     nelts = 0
#     @inbounds for i in 1:nargs
#         n = length(A[i])
#         S[i] = nelts
#         L[i] = n
#         nelts += n
#     end
#     S = KnetArray(S)
#     L = KnetArray(L)
#     X = KnetArray([map(pointer,A)...])
#     Y = KnetArray{$T}(nelts)
#     @knet8($F,(Cint,Ptr{Cint},Ptr{Cint},Ptr{Ptr{$T}},Ptr{$T}),nargs,S,L,X,Y)
#     return Y
# end
# end

# a = [ rand(Float32,200) for i=1:25 ]
# k = map(KnetArray,a)
# @show vcat(a...) == vcat(k...)
# @show vcat(a...) == concat(k...)
# @show @benchmark vcat(a...)
# @show @benchmark vcat(k...)
# @show @benchmark concat(k...)
# @show @benchmark KnetArray(vcat(a...))


text:
#! In this folder I collect functions that are for the version of the model with depreciating information.

"""
    Calculate prior uncertainty after n periods of a stable data regulation regime.
"""
function CalcPriorUncertainty(n, past_posterior; par::ParsGrowth)
    if n == 0
        par.σ²_θ
    else
        par.ρ^2 * past_posterior + par.σ²_η
    end
end

function CalcPosteriorUncertainty(past_prior, par::ParsGrowth)
    #1/(past_prior^-1 + par.σ²_ϵ^-1 + par.σ²_η^-1)
    1/(past_prior^-1 + par.σ²_ϵ^-1)
end

function CalcAllUncertainty(n_max; par::ParsGrowth)
    prior_unc = zeros(n_max + 1)
    post_unc = zeros(n_max + 1)
    for n in 0:n_max
        if n == 0
            prior_unc[n+1] = CalcPriorUncertainty(n, 0; par)
        else
            prior_unc[n+1] = CalcPriorUncertainty(n, post_unc[n]; par)
        end
        post_unc[n+1] = CalcPosteriorUncertainty(prior_unc[n+1], par::ParsGrowth)
    end

    return prior_unc, post_unc
end

function Transition_Probability_dep(par::ParsGrowth, bpar::BParsGrowth; normalize::Symbol = :yes)
    @unpack_ParsGrowth par
    @unpack_BParsGrowth bpar
    Π    = zeros(N, N)

    # Rows: Today's state i
    # Columns: Tomorrow's state j
    for i = 1:N, j = 1:N
        if j == 1
            d = a_grid[j+1] - a_grid[j]

            c_plus = (σ²_ϵ * ((V_n^-1 + σ²_ϵ^-1) * (a_grid[j] + d / 2) / ρ - V_n^-1 * a_grid[i]) - a_grid[i]) / √(V_n + σ²_ϵ)

            Π[i, j] = normcdf(c_plus)
        elseif j == N
            d = a_grid[j] - a_grid[j-1]

            c_minus = (σ²_ϵ * ((V_n^-1 + σ²_ϵ^-1) * (a_grid[j] - d / 2) / ρ - V_n^-1 * a_grid[i]) - a_grid[i]) / √(V_n + σ²_ϵ)

            Π[i, j] = 1 - normcdf(c_minus)
        else
            d_up = a_grid[j+1] - a_grid[j]
            d_down = a_grid[j] - a_grid[j-1]

            c_plus = (σ²_ϵ * ((V_n^-1 + σ²_ϵ^-1) * (a_grid[j] + d_up / 2) / ρ - V_n^-1 * a_grid[i]) - a_grid[i]) / √(V_n + σ²_ϵ)

            c_minus = (σ²_ϵ * ((V_n^-1 + σ²_ϵ^-1) * (a_grid[j] - d_down / 2) / ρ - V_n^-1 * a_grid[i]) - a_grid[i]) / √(V_n + σ²_ϵ)

            Π[i, j] = normcdf(c_plus) - normcdf(c_minus)
        end
    end
    # Renormalize the transition matrix
    if normalize == :yes
        for i = 1:N
            Π[i, :] = Π[i, :] ./ sum(Π[i, :])
        end
    end
    return Π
end


text:
# Try to approximate squared-exponential kernel:
using PyPlot
using Optim
using Roots
using Polynomials
include("sturms_theorem.jl")

#function gaussian_approx6()

ntime = 1000
time = collect(linspace(0,10,ntime))
gaussian=exp(-time.^2/2.)

function gauss_model(x)
model = zeros(ntime)
for k=1:nexp
  if k > 1
    x[2+(k-1)*4]=0.0
    x[4+(k-1)*4]=0.0
  end
  model += x[1+(k-1)*4].*exp(-x[3+(k-1)*4].*time).*cos(x[4+(k-1)*4].*time)
  model += x[2+(k-1)*4].*exp(-x[3+(k-1)*4].*time).*sin(x[4+(k-1)*4].*time)
end
return model
end

function compute_chi(x)
model = gauss_model(x)
sum_coeff = zeros(x)
# Compute Taylor-series coefficients:
for k=1:nexp
  sum_coeff[1] += x[1+(k-1)*4]
  sum_coeff[2] +=-x[1+(k-1)*4]*x[3+(k-1)*4]+x[2+(k-1)*4]*x[4+(k-1)*4]
  sum_coeff[3] += x[1+(k-1)*4]*x[3+(k-1)*4]^2-2*x[2+(k-1)*4]*x[3+(k-1)*4]*x[4+(k-1)*4]-x[1+(k-1)*4]*x[4+(k-1)*4]^2
  sum_coeff[4] +=-x[1+(k-1)*4]*x[3+(k-1)*4]^3+3*x[2+(k-1)*4]*x[3+(k-1)*4]^2*x[4+(k-1)*4]+3*x[1+(k-1)*4]*x[3+(k-1)*4]*x[4+(k-1)*4]^2-x[2+(k-1)*4]*x[4+(k-1)*4]^3
end

return sum((gaussian-model).^2) +100.*((sum_coeff[1]-1.0)^2+sum_coeff[2]^2+(sum_coeff[3]+1.0)^2+sum_coeff[4]^2)+(sturms_theorem(x) > 0)*1e10
end

nexp = 4
npar = 4*nexp

#atrial = [-0.5185744409672999,0.2854238724264553,1.3848025999389613,1.8339730386273703,63.11883856322811,-0.002585613714221086,1.479593361171615,-0.001569459359505545,-32.174013349165506,-0.00439276758185438,1.593669981824972,0.00035151389750538633,-29.42595216593894,-0.0013496416338785663,1.3889097102817822,0.002608575323797371]

#atrial = [-0.5185618444470842,0.28542350796810384,1.3848024403617407,1.8339737219232786,63.118858757909,-0.002595606253169336,1.4795927335503265,-0.001568656597234968,-32.17399733592503,-0.004866618683786094,1.5936706619228656,0.0015491916728168549,-29.425929677155498,-0.0013496810594145385,1.38891017342687,0.002608476149940659,0.0006171992813012984,0.001000720300704954,0.10003813140421611,0.099854716608336]
atrial = [-0.5185618444470842,0.28542350796810384,1.3848024403617407,1.8339737219232786,63.118858757909,-0.002595606253169336,1.4795927335503265,-0.001568656597234968,-32.17399733592503,-0.004866618683786094,1.5936706619228656,0.0015491916728168549,-29.425929677155498,-0.0013496810594145385,1.38891017342687,0.002608476149940659]

#chibest=sum((gaussian-model).^2)
chibest=compute_chi(atrial)

println("Initial chi-square: ",chibest)
read(STDIN,Char)

#result = optimize(compute_chi, atrial, BFGS(), OptimizationOptions(autodiff = true))
result = optimize(compute_chi, atrial)

println(Optim.minimizer(result),Optim.minimum(result))
abest= Optim.minimizer(result)
clf()
model = gauss_model(abest)
plot(time,gaussian)
plot(time,model)
plot(time,(gaussian-model))
chi=sum((gaussian-model).^2)
println("New minimum: ",chi," ",abest," ",std(gaussian-model)," ",maximum(abs(gaussian-model)))
read(STDIN,Char)
# Now, plot the power spectrum:
nfreq = 100000
#omega = logspace(-1,1,nfreq)
omega = logspace(-3,1,nfreq)
psd_grid = zeros(nfreq)
clf()
#plot(omega,psd_grid)
# Compute the asymptotic value:
for k=1:nexp
  a1 = abest[1+(k-1)*4]
  a2 = abest[2+(k-1)*4]
  b  = abest[3+(k-1)*4]
  om = abest[4+(k-1)*4]
  psdt = (a2*om.*(b^2 + om^2 - omega.^2) + a1*b.*(b^2 + om^2 + omega.^2))./
         (b^4 + (om^2 - omega.^2).^2 + 2*b^2.*(om^2 + omega.^2))
  psd_grid += psdt
end
clf()
plot(omega,psd_grid)
plot(omega,exp(-omega.^2./2)*sqrt(pi/2))
plot(omega,psd_grid-exp(-omega.^2./2)*sqrt(pi/2))
#return
#end

max_stars_repo_name: dfm/ess

max_issues_repo_name: dfm/ess

max_forks_repo_name: dfm/ess


text:
using ControlSystems
using Statistics

"""
    sv, sx, sigma_se, max_sx, max_sy = simulateMS(v, M=16, mtf=zpk(1,0,1,1), d=0, dw=[1..], sx0=[0..])

Simulate the vector-based Mismatch-Shaping for a multi-element DAC.
Inputs:
- `v`:        A vector of the digital input values. `v` is in (-`M`:2:`M`) if `dw`=[1..]
              otherwise `v` is in [-sum(`dw`),sum(`dw`)].
- `M`:        The number of elements.
- `mtf`:	  The mismatch-shaping transfer function, given in zero-pole form.
- `d`:        Dither uniformly distributed in [-`d`,`d`] is added to sy.
- `dw`:       A vector of dac element weights
- `sx0`:      A matrix whose columns are the initial states of the ESL.
Outputs:
- `sv`:       An MxN matrix whose columns are the selection vectors.
- `sx`:       An orderxM matrix containing the final state of the ESL.
- `sigma_se`: The rms value of the selection error.
- `max_sx`:   The maximum absolute value of the state for all modulators.
- `max_sy`:   The maximum absolute value of the input to the VQ.

For a brief description of the theory of mismatch-shaping DACs, see
R. Schreier and B. Zhang "Noise-shaped multibit D/A convertor employing
unit elements," Electronics Letters, vol. 31, no. 20, pp. 1712-1713,
Sept. 28 1995.
"""
function simulateMS(v, M=16, mtf=zpk([1], [0], 1, 1), d=0, dw=[], sx0=[])
    order = length(mtf.matrix[1].p)
    if isempty(sx0)
        sx0 = zeros(order, M)
    end
    if isempty(dw)
        dw = ones(M)
    end

    # B/A = MTF-1
    mtf = tf(mtf)
    num = mtf.matrix[1].num[end:-1:0]
    den = mtf.matrix[1].den[end:-1:0]
    A = real(-den[2:order+1])
    B = real(num[2:order+1] + A)

    N = length(v)
    sv = zeros(Int, M, N)

    sx = sx0'
    max_sx = maximum(abs.(sx))
    max_sy = 0
    sum_se2 = 0

    for i = 1:N
        # Compute the sy vector.
        sy = sx * B
        # Normalize sy for a minimum value of zero.
        sy = sy .- minimum(sy)
        # Pick the elements that have the largest desired_usage (sy) values.
        sv[:,i] = selectElement(v[i], sy + d * (2*rand(M) .- 1), dw)
        # Compute the selection error.
        se = sv[:,i] - sy
        # Compute the new sx matrix
        sxn = sx * A + se
        sx = [sxn sx[:,1:order-1]]
        # Keep track of some statistics.
        sum_se2 = sum_se2 + sum((se .- mean(se)).^2)
        max_sx = max(max_sx, maximum(abs.(sxn)))
        max_sy = max(max_sy, maximum(abs.(sy)))
    end
    sigma_se = sqrt(sum_se2 / (M*N))

    return sv, sx', sigma_se, max_sx, max_sy
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
export MaximalIndependentSet

struct MaximalIndependentSet end

"""
    independent_set(g, MaximalIndependentSet())

Find a random set of vertices that are independent (no two vertices are adjacent to each other) and
it is not possible to insert a vertex into the set without sacrificing the independence property.

### Implementation Notes
Performs [Approximate Maximum Independent Set](https://en.wikipedia.org/wiki/Maximal_independent_set#Sequential_algorithm) once.
Returns a vector of vertices representing the vertices in the independent set.

### Performance
Runtime: O(|V|+|E|)
Memory: O(|V|)
Approximation Factor: maximum(degree(g))+1
"""
function independent_set(
    g::AbstractGraph{T},
    alg::MaximalIndependentSet
    ) where T <: Integer

    ind_set = Vector{T}()
    sizehint!(ind_set, nv(g))
    deleted = falses(nv(g))

    for v in shuffle(vertices(g))
        deleted[v] && continue

        deleted[v] = true
        push!(ind_set, v)
        deleted[neighbors(g, v)] .= true
    end

    return ind_set
end


text:
"""
    Base.convert(::Type{Matrix{Char}}, strings::Vector{String})

Converts a `Vector{String}` to a `Matrix{Char}`, where each column contains
the characters from a string in `strings`, truncated to the length of the
shortest string.
"""
function Base.convert(::Type{Matrix{Char}}, strings::Vector{String})
    sorted = sort!(strings, by = length)
    max_length = length(sorted[1])
    return reduce(hcat, collect(s)[1:max_length] for s in strings)
end

"Helper function to identify iterables where all the values are the same"
all_equal(v) = all(x -> x == v[1], v)

"""
    longest_common_prefix(strings::Vector{String})

Given a list of strings, returns a string containing the longest common prefix
among all strings. Returns an empty string if no matching prefix is found.
"""
function longest_common_prefix(strings::Vector{String})
    char_matrix    = convert(Matrix{Char}, strings)
    matching_chars = Iterators.takewhile(all_equal, eachrow(char_matrix))
    return map(x -> first(x), matching_chars) |> join
end

# Tests! -----------------------------------------------------------------------

using Test

@testset "Should return any single letter as-is" begin
    @test longest_common_prefix(["a"]) == "a"
    @test longest_common_prefix(["b"]) == "b"
    @test longest_common_prefix(["z"]) == "z"
end

@testset "Should identify single-letter prefixes" begin
    input = ["apple", "aardvark"]
    @test longest_common_prefix(input) == "a"

    input = ["garbage", "graendal", "gift"]
    @test longest_common_prefix(input) == "g"

    input = ["eagle", "easter", "easy", "egg"]
    @test longest_common_prefix(input) == "e"
end

@testset "Should identify multiple-letter prefixes" begin
    input = ["apple", "apricot"]
    @test longest_common_prefix(input) == "ap"

    input = ["garbage", "garden", "gargantuan"]
    @test longest_common_prefix(input) == "gar"

    input = ["castle", "casting", "castanet", "cast"]
    @test longest_common_prefix(input) == "cast"
end

@testset "Should identify when there is no common prefix" begin
    input = ["apple", "banana"]
    @test longest_common_prefix(input) == ""

    input = ["garbage", "garden", "helicopter"]
    @test longest_common_prefix(input) == ""

    input = ["castle", "antimony", "castanet", "antiquity"]
    @test longest_common_prefix(input) == ""
end

@testset "Should identify repeated instances of a single word" begin
    input = ["apple", "apple"]
    @test longest_common_prefix(input) == "apple"

    input = ["apple", "apple", "apple"]
    @test longest_common_prefix(input) == "apple"
end

@testset "Should work on large inputs" begin
    input = fill("variadic", 10000)
    @test longest_common_prefix(input) == "variadic"

    append!(input, fill("variable", 10000))
    append!(input, fill("variety", 10000))
    @test longest_common_prefix(input) == "vari"

    push!(input, "a")
    @test longest_common_prefix(input) == ""
end


text:
# Xpress model attributes

############################################
#
#   Low level attribute getters
#
############################################
"""
    get_intattr(model::Model, ipar::Integer)

Return integer value corresponding to attribute with number `ipar`
"""
function get_intattr(model::Model, ipar::Integer)
    out = Array{Cint}(undef, 1)

    ret = @xprs_ccall(getintattrib, Cint,
                      (Ptr{Nothing}, Cint, Ptr{Cint}),
                      model.ptr_model, Cint(ipar), out)
    if ret != 0
        throw(XpressError(model, ret))
    end

    return convert(Cint, out[1])
end

"""
    get_dblattr(model::Model, ipar::Integer)

Return double value corresponding to attribute with number `ipar`
"""
function get_dblattr(model::Model, ipar::Integer)
    out = Array{Float64}(undef, 1)

    ret = @xprs_ccall(getdblattrib, Cint,
                      (Ptr{Nothing}, Cint, Ptr{Float64}),
                      model.ptr_model, Cint(ipar), out)
    if ret != 0
        throw(XpressError(model, convert(Int, ret)))
    end

    return out[1]::Float64
end

"""
    get_strattr(model::Model, ipar::Integer)

Return string value corresponding to attribute with number `ipar`
"""
function get_strattr(model::Model, ipar::Integer)
    out = zeros(Cchar,256)

    ret = @xprs_ccall(getstrattrib, Cint,
                      (Ptr{Nothing}, Cint, Ptr{Cchar}),
                      model.ptr_model, Cint(ipar), out)
    if ret != 0
        throw(XpressError(model, convert(Int, ret)))
    end

    return unsafe_string(pointer(out))
end

############################################
#
#   Macros for array definition
#
############################################

macro xprs_int_attr(fun, attrname)
    @eval $(fun)(model::Model) = get_intattr(model, $(attrname))
end

macro xprs_dbl_attr(fun, attrname)
    @eval $(fun)(model::Model) = get_dblattr(model, $(attrname))
end

macro xprs_str_attr(fun, attrname)
    @eval $(fun)(model::Model) = get_strattr(model, $(attrname))
end

############################################
#
#   Model attributes
#
############################################

# basic attributes

@xprs_int_attr num_presolvedvars     XPRS_COLS
@xprs_int_attr num_presolvedconstrs  XPRS_ROWS
@xprs_int_attr num_sos      XPRS_SETS
@xprs_int_attr num_qconstrs XPRS_QCONSTRAINTS
@xprs_int_attr num_cnzs     XPRS_ELEMS
@xprs_int_attr num_qnzs     XPRS_QELEMS
@xprs_int_attr num_qcnzs    XPRS_QCELEMS

@xprs_int_attr num_intents  XPRS_MIPENTS

@xprs_int_attr num_setmembers XPRS_SETMEMBERS

@xprs_dbl_attr obj_sense    XPRS_OBJSENSE

@xprs_int_attr num_originalvars     XPRS_ORIGINALCOLS
@xprs_int_attr num_vars     XPRS_ORIGINALCOLS
@xprs_int_attr num_originalconstrs  XPRS_ORIGINALROWS
@xprs_int_attr num_constrs  XPRS_ORIGINALROWS

# derived attribute functions

"""
    num_linconstrs(model::Model)

Return the number of purely linear contraints in the Model
"""
num_linconstrs(model::Model) = num_constrs(model) - num_qconstrs(model)

"""
    model_sense(model::Model)

Return a symbol that encodes the objective function sense.
The output is either `:minimize` or `:maximize`
"""
model_sense(model::Model) = obj_sense(model) == XPRS_OBJ_MINIMIZE ? (:minimize) : (:maximize)

"""
    is_qcp(model::Model)

Return `true` if there are quadratic  constraints in the Model
"""
is_qcp(model::Model) = num_qconstrs(model) > 0

"""
    is_mip(model::Model)

Return `true` if there are integer entities in the Model
"""
is_mip(model::Model) = (num_intents(model)+num_sos(model)) > 0

"""
    is_qp(model::Model)

Return `true` if there are quadratic terms inthe objective in the Model
"""
is_qp(model::Model) = num_qnzs(model) > 0

"""
    model_type(model::Model)

Return a symbol enconding the type of the model.]
Options are: `:LP`, `:QP` and `:QCP`
"""
function model_type(model::Model)
    is_qp(model)  ? (:QP)  :
    is_qcp(model) ? (:QCP) : (:LP)
end

############################################
#
#   The show method for model
#
#   - Based on attributes
#
############################################

"""
    show(io::IO, model::Model)

Prints a simplified model description
"""
function Base.show(io::IO, model::Model)
    if model.ptr_model != C_NULL
        println(io, "Xpress Model:"     )# $(model_name(model))")
        if is_mip(model)
            println(io, "    type   : $(model_type(model)) (MIP)")
        else
            println(io, "    type   : $(model_type(model))")
        end
        println(io, "    sense  : $(model_sense(model))")
        println(io, "    number of variables             = $(num_vars(model))")
        println(io, "    number of linear constraints    = $(num_linconstrs(model))")
        println(io, "    number of quadratic constraints = $(num_qconstrs(model))")
        println(io, "    number of sos constraints       = $(num_sos(model))")
        println(io, "    number of non-zero coeffs       = $(num_cnzs(model))")
        println(io, "    number of non-zero qp objective terms  = $(num_qnzs(model))")
        println(io, "    number of non-zero qp constraint terms = $(num_qcnzs(model))")
        println(io, "    number of integer entities = $(num_intents(model))")
    else
        println(io, "Xpress Model: NULL")
    end
end

############################################
#
#   more setters and getters
#
############################################
"""
    set_sense!(model::Model, sense::Symbol)

Set the sense of the model.
Options are `:minimize` and `:maximize`
"""
function set_sense!(model::Model, sense::Symbol)
    v = sense == :maximize || sense == :Max ? XPRS_OBJ_MAXIMIZE :
        sense == :minimize || sense == :Min ? XPRS_OBJ_MINIMIZE :
        throw(ArgumentError("Invalid model sense."))

    set_sense!(model, v)

    return nothing
end
function set_sense!(model::Model, sense::Integer)

    ret = @xprs_ccall(chgobjsense, Cint, (
            Ptr{Nothing},    # model
            Cint          # sense
            ),
            model.ptr_model, Cint(sense))

    if ret != 0
        throw(XpressError(model))
    end

    return nothing
end

"""
    set_objcoeffs!{I<:Integer,R<:Real}(model::Model, inds::Vector{I}, c::Vector{R})

Sets coefficients `c` given indices in `inds` in the objective function of `model`
"""
function set_objcoeffs!(model::Model, inds::Vector{I}, c::Vector{R}) where {I<:Integer,R<:Real}
    # n = num_vars(model)
    # _chklen(c,n)
    _chklen(inds, length(c))

    ret = @xprs_ccall(chgobj, Cint, (
        Ptr{Nothing},     # model
        Cint,          # nels
        Ptr{Cint}, # inds
        Ptr{Float64} # vals
        ),
        model.ptr_model, Cint(length(c)), ivec(inds) .- Cint(1), fvec(c))

    if ret != 0
        throw(XpressError(model))
    end

    return nothing
end
"""
    set_objcoeffs!{R<:Real}(model::Model, c::Vector{R})
    set_obj!(model, c)

Sets coefficients in the objective of the dirst `length(c)` all variables.
"""
set_objcoeffs!(model::Model, c::Vector{R}) where {R<:Real} = set_objcoeffs!(model, collect(1:length(c)), c)
set_obj!(model, c) = set_objcoeffs!(model,c)

"""
    get_lb!(model::Model, lb::Vector{Float64})

Return the lower bounds for all variables in the vector lb.
"""
function get_lb!(model::Model, lb::Vector{Float64}, colb::Integer, cole::Integer)

    _chklen(lb,cole-colb+1)

    ret = @xprs_ccall(getlb, Cint, (
        Ptr{Nothing},    # model
        Ptr{Float64},
        Cint,
        Cint
        ),
        model.ptr_model, lb, Cint(colb-1), Cint(cole-1))

    if ret != 0
        throw(XpressError(model))
    end

    return nothing
end
function get_lb!(model::Model, lb::Vector{Float64})

    cols = num_vars(model)
    get_lb!(model, lb, 1, cols)

    return nothing
end

"""
    get_lb(model::Model)
    lowerbounds(model::Model)

Return vector of lowebounds with length equals to the number of variables in the model.
"""
function get_lb(model::Model, colb::Integer, cole::Integer)

    out = Array{Float64}(undef, cole-colb+1)

    get_lb!(model, out, colb, cole)

    return out
end
function get_lb(model::Model)

    cols = num_vars(model)
    out = Array{Float64}(undef, cols)

    get_lb!(model, out)

    return out
end
lowerbounds(model::Model) = get_lb(model)

"""
    get_ub!(model::Model, out::Vector{Float64})

Return the upper bounds for all variables in the vector out.
"""
function get_ub!(model::Model, out::Vector{Float64}, colb::Integer, cole::Integer)

    _chklen(out, cole-colb+1)

    ret = @xprs_ccall(getub, Cint, (
        Ptr{Nothing},    # model
        Ptr{Float64},
        Cint,
        Cint
        ),
        model.ptr_model, out, Cint(colb-1), Cint(cole-1))

    if ret != 0
        throw(XpressError(model))
    end

    return nothing
end
function get_ub!(model::Model, ub::Vector{Float64})

    cols = num_vars(model)
    get_ub!(model, ub, 1, cols)

    return nothing
end

"""
    get_ub(model::Model)
    upperbounds(model::Model)
    get_ub(model::Model, colb::Integer, cole::Integer)

Return vector of upperbounds with length equals to the number of variables in the model.
"""
function get_ub(model::Model, colb::Integer, cole::Integer)

    out = Array{Float64}(undef, cole-colb+1)

    get_ub!(model, out, colb, cole)

    return out
end
function get_ub(model::Model)

    cols = num_vars(model)
    out = Array{Float64}(undef, cols)

    get_ub!(model, out)

    return out
end
upperbounds(model::Model) = get_ub(model)

"""
    get_obj!(model::Model, obj::Vector{Float64})

Return the upper bounds for all variables in the vector obj.
"""
function get_obj!(model::Model, obj::Vector{Float64})

    cols = num_vars(model)
    _chklen(obj,cols)

    ret = @xprs_ccall(getobj, Cint, (
        Ptr{Nothing},    # model
        Ptr{Float64},
        Cint,
        Cint
        ),
        model.ptr_model, obj, Cint(0), Cint(cols - 1))

    if ret != 0
        throw(XpressError(model))
    end

    return nothing
end

"""
    get_obj(model::Model)
    upperbounds(model::Model)

Return vector of objective coefficients with length equals to the number of variables in the model.
"""
function get_obj(model::Model)

    cols = num_vars(model)
    out = Array{Float64}(undef, cols)

    get_obj!(model, out)

    return out
end
objcoeffs(model::Model) = get_obj(model)

"""
    get_rhs!(model::Model, rhs::Vector{Float64})

Return the rhs for all constraints in the vector obj.
"""
function get_rhs!(model::Model, out::Vector{Float64}, rowb::Integer, rowe::Integer)

    _chklen(out, rowe-rowb+1)

    ret = @xprs_ccall(getrhs, Cint, (
        Ptr{Nothing},    # model
        Ptr{Float64},
        Cint,
        Cint
        ),
        model.ptr_model, out, Cint(rowb-1), Cint(rowe-1))

    if ret != 0
        throw(XpressError(model))
    end

    return nothing
end
function get_rhs!(model::Model, out::Vector{Float64})

    rows = num_constrs(model)
    get_rhs!(model, out, 1, rows)

    return nothing
end

"""
    get_rhs(model::Model)

Return a vector of rhs with length equals to the number of variables in the model.
"""
function get_rhs(model::Model, rowb::Integer, rowe::Integer)

    out = Array{Float64}(undef, rowe-rowb+1)

    get_rhs!(model, out, rowb, rowe)

    return out
end
function get_rhs(model::Model)

    rows = num_constrs(model)
    out = Array{Float64}(undef, rows)

    get_rhs!(model, out)

    return out
end

"""
    get_rowtype!(model::Model, sense::Vector{Cchar})

Return the sense for all constraints in the vector sense of Cchar.
Options are XPRS_LEQ = 'L', XPRS_GEQ = 'G', XPRS_EQ = 'E'.
"""
function get_rowtype!(model::Model, sense::Vector{Cchar})

    rows = num_constrs(model)
    _chklen(sense,rows)

    ret = @xprs_ccall(getrowtype, Cint, (
        Ptr{Nothing},    # model
        Ptr{Cchar},
        Cint,
        Cint
        ),
        model.ptr_model, sense, Cint(0), Cint(rows.-1))

    if ret != 0
        throw(XpressError(model))
    end

    return nothing
end

"""
    get_rowtype(model::Model)

Return a vector of senses with length equals to the number of constraints in the model.
Options are XPRS_LEQ = 'L', XPRS_GEQ = 'G', XPRS_EQ = 'E'.
"""
function get_rowtype(model::Model)

    rows = num_constrs(model)
    out = Array{Cchar}(undef, rows)

    get_rowtype!(model, out)

    return out
end

const sense_map = Dict{Cchar, Symbol}(
    XPRS_LEQ => :<=,
    XPRS_GEQ => :>=,
    XPRS_EQ => :(==),
)

get_sense(model::Model) = map(x->sense_map[x], get_rowtype(model))

"""
    get_coltype!(model::Model, coltype::Vector{Cchar})

Return the type for all constraints in the vector coltype of Cchar.
Options are XPRS_CONTINUOUS = 'C', XPRS_INTEGER = 'I', XPRS_BINARY = 'B'.
"""
function get_coltype!(model::Model, coltype::Vector{Cchar})
    cols = num_vars(model)
    _chklen(coltype,cols)

    ret = @xprs_ccall(getcoltype, Cint, (
        Ptr{Nothing},    # model
        Ptr{Cchar},
        Cint,
        Cint
        ),
        model.ptr_model, coltype, 0, cols - Cint(1))

    if ret != 0
        throw(XpressError(model))
    end

    return nothing
end

"""
    get_coltype(model::Model)

Return a vector of type with length equals to the number of variables in the model.
Options are XPRS_CONTINUOUS = 'C', XPRS_INTEGER = 'I', XPRS_BINARY = 'B'.
"""
function get_coltype(model::Model)

    cols = num_vars(model)
    out = Array{Cchar}(undef, cols)

    get_coltype!(model, out)

    return out
end

function unsafe_chgbounds!(model::Model, len::Cint, inds::Vector{Cint}, btype::Vector{Cchar},  lb::Vector{Float64})

    ret = @xprs_ccall(chgbounds, Cint, (
        Ptr{Nothing},    # model
        Cint,
        Ptr{Cint},
        Ptr{Cchar},
        Ptr{Float64}
        ),
        model.ptr_model, len, inds, btype, lb)

    if ret != 0
        throw(XpressError(model))
    end
    return nothing
end
chgbounds!(model::Model, inds::Vector{Cint}, btype::Vector{Cchar},  lb::Vector{Float64}) = unsafe_chgbounds!(model, Cint(length(inds)), inds .- Cint(1), btype,  lb)

"""
    set_lb!{I<:Integer, R<:Real}(model::Model, inds::Vector{I}, lb::Vector{R})

Sets lower bounds `lb` given variable indices in `inds` of `model`

    set_lb!{R<:Real}(model::Model, lb::Vector{R})

Sets lower bounds to all variables up to `length(lb)`.
`length(lb)` must be smaller than the number of variables.
"""
function set_lb!(model::Model, inds::Vector{I}, lb::Vector{R}) where {I<:Integer, R<:Real}
    nbnds = length(inds)
    chgbounds!(model, ivec(inds), cvecx('L', nbnds),  fvec(lb))
    return nothing
end
function set_lb!(model::Model, lb::Vector{R}) where R<:Real
    cols = num_vars(model)
    _chklen(lb,cols)
    ind = inds32(cols)
    set_lb!(model, ind, lb)
    return nothing
end

"""
    set_ub!{I<:Integer, R<:Real}(model::Model, inds::Vector{I}, ub::Vector{R})

Sets upper bounds `ub` given variable indices in `inds` of `model`

    set_ub!{R<:Real}(model::Model, ub::Vector{R})

Sets upper bounds to all variables up to `length(ub)`.
`length(ub)` must be smaller than the number of variables.
"""
function set_ub!(model::Model, inds::Vector{I}, ub::Vector{R}) where {I<:Integer, R<:Real}
    nbnds = length(inds)
    chgbounds!(model, ivec(inds), cvecx('U', nbnds),  fvec(ub))
    return nothing
end
function set_ub!(model::Model, ub::Vector{R}) where R<:Real
    cols = num_vars(model)
    _chklen(ub,cols)
    ind = inds32(cols)
    set_ub!(model,ind,ub)
    return nothing
end

function unsafe_chgrhs!(model::Model, nels::Cint, inds::Vector{Cint}, rhs::Vector{Float64})
    ret = @xprs_ccall(chgrhs, Cint, (
        Ptr{Nothing},    # model
        Cint,
        Ptr{Cint},
        Ptr{Float64}
        ),
        model.ptr_model, nels, inds, rhs)

    if ret != 0
        throw(XpressError(model))
    end
    return nothing
end
chgrhs!(model::Model, inds::Vector{Cint}, rhs::Vector{Float64}) = unsafe_chgrhs!(model, Cint(length(inds)), inds .- Cint(1), rhs)

"""
    set_rhs!{I<:Integer, R<:Real}(model::Model, inds::Vector{I}, rhs::Vector{R})

Sets coefficients `rhs` given indices in `inds` in the rhs of `model`

    set_rhs!{R<:Real}(model::Model, lb::Vector{R})

Sets coefficients in the rhs of all constraints up to `length(rhs)`.
`length(rhs)` must be smaller than the number of constraints.
"""
function set_rhs!(model::Model, inds::Vector{I}, rhs::Vector{R}) where {I<:Integer, R<:Real}
    nels = length(inds)
    chgrhs!(model, ivec(inds), fvec(rhs))
    return nothing
end
function set_rhs!(model::Model, rhs::Vector{R}) where R<:Real
    rows = length(rhs)
    set_rhs!(model, inds32(rows), rhs)
end

"""
    set_rowtype!{I<:Integer}(model::Model, inds::Vector{I}, senses::Vector{Cchar})

Sets row type `senses` for given indices in `inds`

    set_rowtype!(model::Model, senses::Vector{Cchar})

Sets row type in all constraints up to `length(senses)`.
`length(senses)` must be smaller than the number of constraints.
"""
function set_rowtype!(model::Model, inds::Vector{I}, senses::Vector{Cchar}) where I<:Integer

    rows = length(senses)
    ret = @xprs_ccall(chgrowtype, Cint, (
        Ptr{Nothing},    # model
        Cint,
        Ptr{Cint},
        Ptr{Cchar}
        ),
        model.ptr_model, Cint(rows) , ivec(inds) .- Cint(1), cvec(senses) )

    if ret != 0
        throw(XpressError(model))
    end

    return nothing
end
function set_rowtype!(model::Model, senses::Vector{Cchar})

    rows = length(senses)
    inds = inds32(rows)
    set_rowtype!(model, inds, senses)

    return nothing
end

"""
    set_constrLB!{R<:Real}(model::Model, lb::Vector{R})

Change constraints lower bounds up to `length(lb)`
"""
function set_constrLB!(model::Model, lb::Vector{R}) where R<:Real
    # should only work for liner constraints

    nlrows = num_linconstrs(model)
    lrows = get_lrows(model)[1:length(lb)]

    senses = get_rowtype(model)
    rhs    = get_rhs(model)
    sense_changed = false

    for i = 1:length(lrows)
        #if senses[lrows[i]] == XPRS_GEQ || senses[lrows[i]] == XPRS_EQ
            # Do nothing
        if senses[lrows[i]] == XPRS_LEQ
            if lb[i] != -Inf
                # LEQ constraint with non-NegInf LB implies a range
                if isapprox(lb[i], rhs[lrows[i]])
                    # seems to be an equality
                    senses[lrows[i]] = XPRS_EQ
                    sense_changed = true
                else
                    error("Tried to set LB != -Inf on a LEQ constraint (index $i)")
                end
            else
                lb[i] = rhs[lrows[i]]
            end
        end
    end
    if sense_changed
        set_rowtype!(model, senses)
    end
    set_rhs!(model, lrows, lb)
end

"""
    set_constrUB!{R<:Real}(model::Model, ub::Vector{R})

Change constraints upper bounds up to `length(ub)`
"""
function set_constrUB!(model::Model, ub::Vector{R}) where R<:Real

    nlrows = num_linconstrs(model)
    lrows = get_lrows(model)[1:length(ub)]

    senses = get_rowtype(model)
    rhs    = get_rhs(model)
    sense_changed = false
    for i = 1:length(lrows)
        #if senses[lrows[i]] == XPRS_LEQ || senses[lrows[i]] == XPRS_EQ
            # Do nothing
        if senses[lrows[i]] == XPRS_GEQ
            if ub[i] != Inf
                # GEQ constraint with non-PosInf UB implies a range
                if isapprox(ub[i], rhs[lrows[i]])
                    # seems to be an equality
                    senses[lrows[i]] = XPRS_EQ
                    sense_changed = true
                else
                    error("Tried to set UB != +Inf on a GEQ constraint (index $i)")
                end
            else
              ub[i] = rhs[lrows[i]]
            end
        end
    end
    if sense_changed
        set_rowtype!(model, senses)
    end
    set_rhs!(model, lrows, ub)
end


text:
const population = 5_837_213

struct Data
    turing_data
    cases
    deaths
    hospit
    num_observations
    num_total
    dates
    pop
    predictors
    seroprev
    lockdown
    observations_end
    cases_start
    startdates
end

dan2eng = Dict(
    "Dato"          => "date",
    "Total"         => "total",
    "Hovedstaden"   => "capital",
    "Sj\xe6lland"   => "zealand",
    "Sjælland"      => "zealand",
    "Syddanmark"    => "south",
    "Midtjylland"   => "central",
    "Nordjylland"   => "north",
    "Ukendt Region" => "unknown"
)

function missing2zero(xs)
    i = findfirst(!ismissing, xs)
    xs[1:i] .= 0
    convert(Vector{Int64}, xs)
end
# weekly_average(xs) = [mean(xs[i:(i+6)]) for i in 1:div(length(xs),7)]
# weekly_average(deaths)

function padlast(xs, l)
    n  = length(xs)
    i  = findlast(x -> typeof(x) <: Real, xs) #avoid catching a missing value
    ys = Vector{Union{Float64, Missing}}(undef, l)
    ys[1:n] = xs
    ys[n+1:end] .= xs[i]
    ys
end

function padmissing(xs, l)
    n  = length(xs)
    ys = Vector{Union{eltype(xs), Missing}}(missing, l)
    ys[1:n] = xs
    ys
end

function cum2counts(x)
    x = x |> missing2zero |> diff
    pushfirst!(x, 0)
    x[x.<0] .= 0 # ignore data inconsistencies
    x
end

const dformat1 = DateFormat("d u Y")
const dformat2 = DateFormat("d-u-Y")
function parsetime(ds)
    ds = replace(ds, "Mrz"=>"Mar")
    ds = replace(ds, "Dez"=>"Dec")
    ds = replace(ds, "."=>"")
    d = try
        Date(ds, dformat1)
    catch
        try
            Date(ds, dformat2)
        catch
            @error "Not a valid time format: $ds"
        end
    end
    d + Year(2000)
end

function znorm!(xs)
    xs .-= mean(skipmissing(xs))
    xs ./= std(skipmissing(xs))
    nothing
end

function minmaxnorm!(xs)
    xs .-= minimum(skipmissing(xs))
    xs ./= maximum(skipmissing(xs))
    nothing
end

function surveynorm!(xs)
    i = findfirst(x->!ismissing(x), xs)
    xs[i] == 0. && return nothing
    xs ./= abs(xs[i])
    xs .-= xs[i]
    return nothing
end

# findall(ismissing, xs[mask])
function preproc(xs::AbstractVector, num_obs, Ntot; normalize=nothing)
    mask = 1:num_obs
    xs = padlast(xs[mask], Ntot)
    !isnothing(normalize) && normalize(xs)
    xs
end

# preproc(dk.Rt, num_obs, Ntot; normalize=surveynorm!)
# function select_predictors(dk, predictors, num_obs, Ntot, surveyvars)
#     predictors = String.(predictors)
#
#     #preprocess mobility data
#     dk[!,:nonresidential_mobility] = preproc(
#         (
#             dk.retail_and_recreation_percent_change_from_baseline +
#             dk.grocery_and_pharmacy_percent_change_from_baseline +
#             dk.workplaces_percent_change_from_baseline
#         ) / 3 / 100 , num_obs, Ntot; normalize=nothing)
#     # begin
#     #     Mavg = (
#     #     preproc(dk.retail_and_recreation_percent_change_from_baseline, num_obs, Ntot) +
#     #     preproc(dk.grocery_and_pharmacy_percent_change_from_baseline, num_obs, Ntot) +
#     #     preproc(dk.workplaces_percent_change_from_baseline, num_obs, Ntot)
#     #     ) ./ -3
#     # end
#     dk[!,:residential_mobility] = preproc(dk.residential_percent_change_from_baseline, num_obs, Ntot; normalize=nothing) ./ 100
#
#     #preprocess survey data
#     for k in surveyvars
#         dk[!,k] = preproc(dk[:,k], num_obs, Ntot; normalize=surveynorm!)
#     end
#
#     dk[!, :lockdown] .= 0.
#     t = findfirst(==(Date("2020-03-14")), dk.date) # lockdown date in DK: 11. (partial), 17. (full)
#     dk[t:end, :lockdown] .+= 1
#     # select predictors
#     miss = .!( in.(predictors, Ref(names(dk)) ) )
#     any(miss) && (@error "predictor $(predictors[miss]) not found. Select either nonresidential_mobility, residential_mobility or from survey: $(surveyvars)")
#     Array( dk[!,predictors] )
# end

function select_predictors(covariates, predictors, num_obs = nothing, num_tot = nothing)
    predictors = String.(predictors)

    cov_names = filter(x->x!=="date", names(covariates))
    for k in cov_names
        if isnothing(num_obs) || isnothing(num_tot)
            surveynorm!(covariates[!,k])
        else
            covariates[!,k] = preproc(covariates[:,k], num_obs, num_tot; normalize=surveynorm!)
        end
    end

    miss = .!( in.(predictors, Ref(cov_names) ) )
    any(miss) && (@error "predictor $(predictors[miss]) not found. Select either nonresidential_mobility, residential_mobility or from survey: $(surveyvars)")
    Array( covariates[!,predictors] )
end

rolling(xs) = begin
    n  = length(xs)
    ys = Vector{Union{Missing, Float64}}(missing, n)
    ys[1:n-6] = [maximum(xs[i:i+6]) for i in 1:n-6]
    ys
end

function gamma(μ, cv)
    α = 1/cv^2
    θ = μ/α
    Gamma(α, θ)
end

function nconvolve(f1, f2, ts)
    nbins = length(ts)
    f(τ, t) = f1(τ) * f2(t - τ)
    f3 = zeros(nbins)
    for (i, t) in enumerate(ts)
        integral, err = quadgk(x->f(x, t), 0, t, rtol=1e-8)
        f3[i] = integral
    end
    return f3
end

function inf2hosp(nbins=50, ihr = 1)
    # SSI: ihr ∈ 0.05%-0.5%
    # derived from Verity et al. See parameters.jl for details
    # incubatine time ~ gamma(5.1, .86) see Flaxman et al
    # onset to hospitalization see Zhang et al. Lancet (2020):
    # "Evolving epidemiology and transmission dynamics of coronavirus disease 2019 outside Hubei province, China: a descriptive and modelling study"
    ts = 1.5:1:(nbins + .5)
    i2c(x) = pdf(gamma(5.1, .86), x)
    # c2h(x) = pdf(gamma(2.6, .92), x) * ihr
    c2h(x) = pdf(gamma(6.49, .92), x) * ihr
    return nconvolve(i2c, c2h, ts)
end

function inf2death(nbins=75, ihr = 1)
    # ifr = 89 / 100_000 for age < 70 see Erikstrup et al. in Clin. Inf. Dis.

    ## report34: onset to death = 19.8, Cv = .85
    ts = 1.5:1:(nbins + .5)
    i2c(x) = pdf(gamma(5.1, .86), x)
    c2d(x) = pdf(gamma(19.8, .85), x) * ihr
    return nconvolve(i2c, c2d, ts)

    # Brauer et al.
    # r = vcat([0], 1.5:1:(nbins + .5))
    # diff(cdf.(NegativeBinomial2(21.82, 14.26), r))

    ## S. Abbott
    # i2o(x) = pdf(LogNormal(1.621, 0.418),x)
    # o2d(x) = pdf(LogNormalMeanStd(13.1, 11.7),x)
    # return nconvolve(i2o, o2d, ts)

    ## Nouvelette
    # GammaMeanStd(18.8, 8.48)
end

function inf2case(nbins = 30)
    ## Imperial Report
    r = vcat([0], 1.5:1:(nbins + .5))
    # return diff(cdf.(gamma(5.1, .86), r))
    ## Brauner et al.
    return diff(cdf.(NegativeBinomial2(10.92, 5.41), r))
    ## S. Abbott
    # i2o(x) = pdf(LogNormal(1.621, 0.418),x)
    # o2c(x) = pdf(LogNormalMeanStd(6.5, 17),x)
    # return nconvolve(i2o, o2c, ts)
end

"""
seropos()[t] is the probablity of being seropositive t days after infection.
For details see imperial report 34
"""
function seropos(nbins=200)
    # λ = 13.3
    # μ = 300 #143.7 - λ
    # κ = 3.67
    # f(x, t, κ, μ, λ) = 1/λ*exp(-x/λ-((t-x)/μ)^κ)
    # c2p(t) = quadgk(x->f(x, t, κ, μ, λ), 0, t, rtol=1e-8)[1]
    # i2c(t) = pdf(gamma(5.1, .86), t)
    # ts = 1.5:1:(nbins + .5)
    # return nconvolve(i2c, c2p, ts)

    i2s = zeros(nbins)
    delay = 14
    i2s[delay:end] .= 1.
    return i2s
end

function serialinterval(nbins = 30)
    # Imperial Report
    r = vcat([0], 1.5:1:(nbins + .5))
    # diff(cdf.(GammaMeanStd(6.5, .62), r))
    # return p / sum(p)
    ## Brauer et al. Science (2020)
    p = diff(cdf.(GammaMeanStd(5.06, 2.11), r))
    return p# / sum(p)
end

# function inf2serocon(nbins)
#     N = 10_000_000
#     # onset to seroconversion time from Imperial Report 34
#     x = rand(gamma(5.1, .86), N) .+ rand(Exponential(18.3), N)
#     h = fit(Histogram, x, vcat([0], 1.5:1:150.5))
#     h.weights / N * ihr
# end

function aggregateregions!(df, mun2reg)
  regions = values(mun2reg) |> unique
  for r in regions
    df[!, r] .= 0
  end
  for m in keys(mun2reg)
    df[!, mun2reg[m]] += df[!,m]
  end
end

# function getcases()
#     municipalities = CSV.File(normpath( homedir(), "data/Statens-Serum-Institut/municipalities.csv"); delim="\t", )|>DataFrame
#     mun2reg = Dict(municipalities.municipality[i] => municipalities.region[i]
#                     for i in 1:size(municipalities,1))
#     cases    = CSV.File(normpath( homedir(), "data/Statens-Serum-Institut/Municipality_cases_time_series.csv") )|>DataFrame
#     cases = select(cases, municipalities.municipality, :date_sample => :date)
#     cases[!, :total] = sum(select(cases, municipalities.municipality)|>eachcol)
#     aggregateregions!(cases, mun2reg)
#     cases
# end

function insert_missing_dates!(df)
    ts = df.date[1] : Day(1) : df.date[end]
    missing_date = .! in.(ts, Ref(df.date))
    if any(missing_date)
        for t in ts[missing_date]
            insertrow!(df, t)
        end
    end
end

function replace_missing!(df, v=0)
    @assert !any( ismissing.( df.date ) )
    for c in eachcol(df)
        replace!(c, missing=>v)
    end
    disallowmissing!(df)
end

function process_ssi(df)
    rename!(df, 1=>:region, 2=>:date, 3=>:counts)
    df.region = replace(df.region, dan2eng...)
    dropmissing!(df, :region)
    df = unstack(df, :date, :region, :counts)
    insert_missing_dates!(df)
    replace_missing!(df)
    DataFrames.transform(df, AsTable(Not(:date)) => sum => :country)
end

# function getdeaths()
#     deaths   = CSV.File(normpath( homedir(), "data/Statens-Serum-Institut/Deaths_over_time.csv"))|>DataFrame
#     deaths = deaths[1:end-1,:]
#     rename!(deaths, "Dato"=>"date", "Antal_døde"=>"total")
#     deaths.date = Date.(deaths.date, Ref(dateformat"y-m-d"))
#     deaths
# end
#
# function gethospit()
#     hosp     = CSV.File(normpath( homedir(), "data/Statens-Serum-Institut/Newly_admitted_over_time.csv"))|>DataFrame
#     dan2eng = Dict(
#         "Dato"          => "date",
#         "Total"         => "total",
#         "Hovedstaden"   => "capital",
#         "Sjælland"      => "zealand",
#         "Syddanmark"    => "south",
#         "Midtjylland"   => "central",
#         "Nordjylland"   => "north",
#         "Ukendt Region" => "unknown"
#     )
#     rename!(hosp, dan2eng)
#     hosp
# end

enddate(df1, df2) = min(df1.date[end], df2.date[end])
enddate(df, args...) = min(df.date[end], enddate(args...))

function padzeros(xs::AbstractVector, padfirst=0, total=length(xs)+padfirst)
    n  = length(xs)
    ys = zeros(eltype(xs), total)
    ys[padfirst+1:n+padfirst] = xs
    ys
end

function insertrow!(df, d::Date, v=0)
    i = findfirst(==(d-Day(1)), df.date)
    v = [eltype(c) == Date ? d : v for c in eachcol(df)]
    insert!.(eachcol(df), i+1, v)
    nothing
end

function padzeros(df::DataFrame, s::Date)
    npad = (df.date[1] - s).value
    dfpad = DataFrame( :date => s:Day(1):df.date[end] )
    for n in names(df)
        n == "date" && continue
        dfpad[!,n] = padzeros(df[!,n], npad)
    end
    dfpad
end

function limit(df::DataFrame, s::Date, e::Date, pad=true)
    pad && df.date[1] > s && ( df = padzeros(df, s) )
    is = findfirst(==(s), df.date)
    ie = findlast(==(e), df.date)
    df[is:ie, :]
end

function parseweek(s)
    _,w = split(s)
    w = parse(Int, w)
    y = w > 40 ? Date("2020") : Date("2021")
    y + Week(w) - Day(3) #center at day 4
end

function parsepercent(x::String)
    x = strip(x, '%')
    x = strip(x)
    x = replace(x, "," => ".")
    parse(Float64, x) / 100
end

function parseCI(x)
    x = strip(x, '%')
    x = replace(x, "," => ".")
    lq, uq = split(x, '-')
    parse(Float64, lq) / 100, parse(Float64, uq) / 100
end

function parsepercent!(df)
    regions = ["capital", "zealand", "south", "central", "north"]
    df.country = parsepercent.(df.country)
    # for r in regions
    #     df[!,r] = parsepercent.(df[:,r])
    # end
end

function findindices(dates, ds)
    is = Int[]
    for d in ds
        i = findfirst(==(d), dates)
        push!(is, i)
    end
    is
end

function readsurvey(fname)
    s = CSV.File(fname)|>DataFrame
    s.date = s.date .|> parsetime
    s.date -= Day(1) # survey setup: contacts in past 24h
    dropmissing!(s)
    for n in names(s)
        c = s[:,n]
        if eltype(c) == String
            c = replace.(c, ","=>".")
            s[!,n] = parse.(Ref(Float64), c)
        end
    end
    s
end
#
# function readsurveys_raw()
#     s1 = normpath( homedir(), "data/covidsurvey/contacts.csv" )  |> readsurvey #projectdir("data", "survey", "timeseries_colleagues_strangers.csv") |> readsurvey
#     s1 = select(s1, Not([:Rt, :uncertainty_lower, :uncertainty_upper]))
#     s2 = normpath( homedir(), "data/covidsurvey/predictors.csv" )  |> readsurvey #projectdir("data", "survey", "predictors.csv") |> readsurvey
#     s = innerjoin(s1, s2, on=:date)
#     s[3:end,:] # skip first for data quality
# end
function readcovariates(; fname=nothing, shift=0, startdate=nothing, enddate=nothing, datecol=:date, kwargs...)
    isnothing(fname) && (@error "please give a file name to covariates_kwargs")
    df = CSV.File(fname)|>DataFrame
    df[!,datecol] += Day(shift)
    !isnothing(startdate) && filter!(row -> row[datecol] >= Date(startdate), df)
    !isnothing(enddate) && filter!(row -> row[datecol] <= Date(enddate), df)
    ts = df[1,datecol] : Day(1) : df[end,datecol]
    @assert all(df[!,datecol] .== ts )
    @assert !any( ismissing.(Array(df)) )
    @assert all( isfinite.(Array(df)) )
    rename!(df, datecol=>:date)
    return df
end

country(df) = DataFrame(:date=>df.date, :country=>df.country)

function holiday(dates)
    specialdays = [Date("2020-12-24"), Date("2020-12-25"), Date("2020-12-31"), Date("2021-01-01")]
    holidays = dates .∈ Ref(specialdays)
    convert(Vector{Int64}, holidays)
end

function observables!(turing_data,
    fname_cases  = normpath( homedir(), "data/Statens-Serum-Institut/dashboard/Regionalt_DB/08_bekraeftede_tilfaelde_pr_dag_pr_regions.csv" ),
    fname_hospit = normpath( homedir(), "data/Statens-Serum-Institut/dashboard/Regionalt_DB/06_nye_indlaeggelser_pr_region_pr_dag.csv" ),
    fname_deaths = normpath( homedir(), "data/Statens-Serum-Institut/dashboard/Regionalt_DB/07_antal_doede_pr_dag_pr_region.csv" )
)

    turing_data["cases"] = CSV.File( fname_cases ) |>
        DataFrame |> process_ssi |> country

    turing_data["hospit"] = CSV.File( fname_hospit ) |>
        DataFrame |> process_ssi |> country

    turing_data["deaths"] = CSV.File( fname_deaths ) |>
        DataFrame |> process_ssi |> country

    return nothing
end

function consistent!(data; epidemicstart = 30)
    @unpack cases, hospit, deaths = data
    data["epidemicstart"] = epidemicstart
    s = data["startdate"]  = startdate(deaths; epidemicstart)
    e = enddate(deaths, cases, hospit)
    data["deaths"] = limit(deaths, s, e)
    data["hospit"] = limit(hospit, s, e)
    data["cases"]  = limit(cases, s, e)
    data["dates"]  = s:Day(1):e
    @assert length(data["dates"]) == size(data["cases"],1) == size(data["hospit"],1) == size(data["deaths"],1)
    return nothing
end

function startdate(df; epidemicstart = 30)
    cumulative = cumsum(df[:,:country])
    start_idx = findfirst(>=(10), cumulative)
    return df.date[start_idx] - Day(epidemicstart)
end

function national_timeseries(dates, startdate, observationsend)
    is = findfirst(==(startdate), dates)
    ie = findfirst(==(observationsend), dates)
    return dates[is:ie]
end

function totaldays(dates, startdates)
    observationsend = last(dates)
    national_dates = national_timeseries(dates, startdates, observationsend)
    return length(national_dates)
end

function df2vec(df, startdate, observationsend)
    vs = Vector{Vector{Int64}}()
    is = findfirst(==(startdate), df.date)
    ie = findfirst(==(observationsend), df.date)
    v = vec(df[is:ie, :country])
    return v
end

function turingformat!(data, observationsend = nothing)
    @unpack dates, cases, hospit, deaths = data
    e = data["observationsend"] = isnothing(observationsend) ? last(dates) : Date(observationsend)
    s = data["startdate"]
    data["num_tot"]       = totaldays(dates, s)
    data["dates_turing"]  = national_timeseries(dates, s, e)
    data["cases_turing"]  = df2vec(cases, s, e)
    data["hospit_turing"] = df2vec(hospit, s, e)
    data["deaths_turing"] = df2vec(deaths, s, e)
    data["num_obs"]       = length(data["cases_turing"])
    return nothing
end

# function cases_start_date!(data, cases_start)
#     @unpack dates_turing, num_obs = data
#     if isnothing(cases_start)
#         data["num_case_obs"] = 0
#         data["cases_start_idx"] = nothing
#     else
#         data["cases_start_idx"] = i = findfirst(==(Date.(cases_start)), dates_turing)
#         data["num_case_obs"] = num_obs - i + 1
#     end
#     return nothing
# end

function covariates!(data, predictors; covariates_kwargs...)
    @unpack dates_turing, num_obs = data
    if isnothing(predictors)
        data["covariates"]       = Array{Float64,2}(undef, 1,1)
        data["num_covariates"]   = 0
        data["covariates_start"] = 0
    else
        covariates       = readcovariates(; covariates_kwargs... )
        covariates       = leftjoin(DataFrame(:date=>dates_turing), covariates, on=:date)
        sort!(covariates, :date)
        covariates_start = findfirst(x->!ismissing(x), covariates[:,Not(:date)][:,1])
        covariates       = covariates[covariates_start:num_obs,:]
        covariates       = select_predictors(covariates, predictors) #num_obs+shift_covariates, num_tot
        covariates       = convert(Array{Float64,2}, covariates)
        data["covariates"]       = covariates
        data["covariates_start"] = covariates_start
        data["num_covariates"]   = size(covariates,2)
    end
end

function lockdown!(data, lockdown)
    @unpack dates_turing = data
    i = findfirst(==(Date(lockdown)), dates_turing)
    data["lockdown_index"] = i
    return nothing
end

function stepindex(n, stepsize)
    index = Vector{Int64}(undef, n)
    step = 1
    for i in 1:n
        index[i] = step
        i % stepsize == 0 && ( step += 1 )
    end
    index
end

function randomwalk!(data, stepsize, covariates_kwargs)
    @unpack num_obs, lockdown_index, covariates_start, num_covariates = data
    n = if covariates_kwargs[:semiparametric] && data["num_covariates"] > 0
        covariates_start-lockdown_index-1
    else
        num_obs - lockdown_index
    end
    i = stepindex(n, stepsize)
    data["num_rt_steps"] = n
    data["rt_step_index"] = i
    return nothing
end

function time_varying_iar!(data, casemodel, iarstep=1)
    @unpack num_obs, num_tot = data
    data["iar_start_idx"] = startindex(data, casemodel)
    data["num_iarsteps"] = num_obs - data["iar_start_idx"] + 1
    # data["iar_idx"] = idx = stepindex(n, stepsize)
    # data["num_iarsteps"] = length(unique(idx[1:num_obs]))
    # iar_idx = get_iar_idx(num_tot, num_obs, cases_start_idx, iarstep)
    # num_iarsteps = length(unique(iar_idx))
    return nothing
end
# ============================================================================
# observation model defaults
@kwdef struct CaseInit <: ObservationInit
    obs_start        ::String  = ""
    obs_stop         ::String  = ""
    delay_dist                 = NegativeBinomial2
    delay_length     ::Int64   = 40
    delay_dispersion ::Float64 = 5.41
    population       ::Int64   = National.population
end

@kwdef struct HospitInit <: ObservationInit
    obs_start        ::String  = ""
    obs_stop         ::String  = ""
    delay_dist                 = NegativeBinomial2
    delay_length     ::Int64   = 40
    delay_dispersion ::Float64 = 5.41
    population       ::Int64   = National.population
end

@kwdef struct DeathInit <: ObservationInit
    obs_start        ::String  = ""
    obs_stop         ::String  = ""
    delay_dist                 = NegativeBinomial2
    delay_length     ::Int64   = 60
    delay_dispersion ::Float64 = 14.26
    population       ::Int64   = National.population
end

# ==============================================================================
# observation model

struct ObservParams{SA,SO,DD,DL,DP,P,H,W}
    start::SA
    stop::SO
    delay_dist::DD
    delay_length::DL
    delay_dispersion::DP
    population::P
    holiday::H
    weekday::W
end

function startindex(data, o)
    dates   = data["dates_turing"]
    if isempty(o.obs_start)
        date = Date(data["startdate"]) + Day(data["epidemicstart"])
        i = findfirst(==(date), dates)
        return i
    else
        date = Date(o.obs_start)
        i = findfirst(==(date), dates)
        return i
    end
end

function stopindex(data, o)
    dates   = data["dates_turing"]
    if isempty(o.obs_stop)
        return data["num_obs"]
    else
        date = Date(o.obs_stop)
        i = findfirst(==(date), dates)
        return i
    end
end

function ObservParams(data, o::ObservationInit)
    params = (
        start            = startindex(data, o),
        stop             = stopindex(data, o),
        delay_dist       = o.delay_dist,
        delay_length     = o.delay_length,
        delay_dispersion = o.delay_dispersion,
        # observations     = observations(data, o),
        population       = o.population,
        holiday          = holiday(data["dates"]),
        weekday          = dayofweek.(data["dates"])
    )
    return ObservParams( params... )
end

ObservParams!(data, o::CaseInit)   = data["casemodel"]   = ObservParams(data, o)
ObservParams!(data, o::HospitInit) = data["hospitmodel"] = ObservParams(data, o)
ObservParams!(data, o::DeathInit)  = data["deathmodel"]  = ObservParams(data, o)
# ObservParams!(data, o::Sero)    = data["seromodel"]   = ObservParams(data, o)

# observations(data, ::CaseInit)   = data["cases_turing"]
# observations(data, ::HospitInit) = data["hospit_turing"]
# observations(data, ::DeathInit)  = data["deaths_turing"]

# ===========================================================================
# sero

# ≈ 15days delay from infection to antibody production
# see Ferguson Report and Zhao 2020: "Antibody responses to SARS-CoV-2 in patients of novel coronavirus disease 2019"
function readsero(;later_dataset = true)
    df0 = CSV.File(projectdir("data/seroprevalence/", "seroprev_early.csv"), delim="\t")|>DataFrame
    df1 = CSV.File(projectdir("data/seroprevalence/", "seroprev.csv"), delim="\t")|>DataFrame
    df1[!,:date] = df1.week .|> parseweek
    select!(df1, Not(:week))
    df = later_dataset ? vcat(df0, df1) : df0
    parsepercent!(df)
    df.CI = df.CI .|> parseCI
    df[!,:std] = [max( abs.(df.country[i] .- df.CI[i])... ) / 1.96 for i in 1:length(df.CI)]
    return df
end

function SeroParams!(data, sero=SeroInit1(); later_dataset = true)
    @unpack dates_turing = data

    df = readsero(;later_dataset)
    s = isempty(sero.obs_start) ? dates_turing[1] : Date(sero.obs_start)
    e = isempty(sero.obs_stop) ? dates_turing[end] : Date(sero.obs_stop)
    df = df[ s .<= df.date .<= e, :]
    data["seroprev"] = df

    # σ = df.std
    # uncertainty_dist = sero.dist.(σ, Ref(sero.cv))
    params = (
        mean  = df.country,
        std   = sero.std,
        CI    = df.CI,
        # dstd  = arraydist( uncertainty_dist),
        index = [findfirst(==(d), dates_turing) for d in df.date],
        dates = df.date,
        delay = sero.delay,
        population = population
    )
    data["sero"] = df.country
    data["seromodel"] = SeroParams(params...)

    return nothing
end

@kwdef struct SeroInit1# <: ObservationInit
    obs_start        ::String  = ""
    obs_stop         ::String  = "2020-07-01"
    delay            ::Int64   = 0
    std              ::Float64 = 0.5
    # population       ::Int64   = National.population
    # dist                       = InverseGamma2
    # cv               ::Float64 = 0.1
end

struct SeroParams{M,S,C,I,A,Y,P}
    mean::M
    std::S
    CIs::C
    # dstd::D
    index::I
    dates::A
    delay::Y
    population::P
end

# function hospitalizations!(data)
#     params = (
#         start            = data["cases_start_idx"],
#         stop             = data["num_observations"]
#         delay_dist       = NegativeBinomial2,
#         delay_num        = 40,
#         delay_dispersion = 5.41,
#         observations     = data["cases_turing"],
#         population       = population,
#     )
#     data["hospit_model"] = Observations( params... )
#     return nothing
# end
#

function load_data(;
    observationsend  = nothing,
    predictors        = nothing,
    rwstep           = 1,
    iarstep          = 1,
    epidemicstart    = 30,
    numimpute        = 6,
    casemodel         = CaseInit(),
    hospitmodel       = HospitInit(),
    deathmodel        = DeathInit(),
    seromodel         = SeroInit1(),
    link              = KLogistic(3.),
    invlink           = KLogit(3.),
    lockdown          = "2020-03-18",
    covariates_kwargs = Dict(
        :semiparametric => false,
        :fname => projectdir("data","smoothed_contact_rates.csv"),
        :shift => -1,
        :startdate => "2020-11-10",
        :enddate => nothing
    )
)

    data = Dict{String, Any}()
    observables!(data)
    consistent!(data; epidemicstart)
    turingformat!(data, observationsend)
    # cases_start_date!(data, cases_start)
    covariates!(data, predictors; covariates_kwargs...)
    lockdown!(data, lockdown)
    randomwalk!(data, rwstep, covariates_kwargs)
    time_varying_iar!(data, casemodel, iarstep)
    ObservParams!(data, casemodel)
    ObservParams!(data, hospitmodel)
    ObservParams!(data, deathmodel)
    SeroParams!(data, seromodel)

    # dk = begin
    #     update && run(`Rscript $(datadir())/COVID19.R`)
    #     dk = load(datadir("COVID19", "DK.rds"))
    #     limit(dk, s, e, false)
    # end

    turing_data = (
        num_impute                  = numimpute,
        num_total_days              = data["num_tot"],
        casemodel                   = data["casemodel"],
        cases                       = data["cases_turing"],
        deathmodel                  = data["deathmodel"],
        deaths                      = data["deaths_turing"],
        epidemic_start              = epidemicstart,
        population                  = population, #denmark statistics 2020 Q4
        serial_interval             = serialinterval(15),#padzeros(serialinterval(30), 0, num_tot), #data_usa.turing_data.serial_intervals[1:50],
        num_si                      = 15,
        lockdown_index              = data["lockdown_index"],
        hospitmodel                 = data["hospitmodel"],
        hospit                      = data["hospit_turing"],
        num_rt_steps                = data["num_rt_steps"],
        rt_step_indices             = data["rt_step_index"],
        iar_start_idx               = data["iar_start_idx"],
        num_iarsteps               = data["num_iarsteps"],
        num_observations            = data["num_obs"],
        link                        = link,
        invlink                     = invlink,
        seromodel                   = data["seromodel"],
        sero                        = data["sero"],
        num_covariates              = data["num_covariates"],
        covariates                  = data["covariates"],
        covariates_start            = data["covariates_start"],
        semiparametric              = covariates_kwargs[:semiparametric],
        rwscale                     = sqrt(rwstep),
        # seroprev_mean               = data["seroprev_mean"],
        # seroprev_std                = data["seroprev_std"],
        # seroprev_idx                = data["seroprev_idx"],

    )

    Data(
        turing_data,
        data["cases"],
        data["deaths"],
        data["hospit"],
        data["num_obs"], #num_observations
        data["num_tot"],
        data["dates"],
        population,
        predictors,
        data["seroprev"],
        lockdown,
        data["observationsend"],
        -1,
        data["startdate"]
    )
end

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
module test_bukdu_head

using Test
using Bukdu
using HTTP: header

struct IndexController <: ApplicationController
    conn::Conn
end

index(c::IndexController) = render(Text, "ok")

routes() do
    get("/", IndexController, index)
end

result = Router.call(Bukdu.head, "/")
@test result.got === "ok"
@test result.resp.status == 200
@test header(result.resp, "Content-Length") == "0"
@test result.resp.body == Vector{UInt8}()

result = Router.call(get, "/")
@test result.got == "ok"
@test result.resp.status == 200
@test header(result.resp, "Content-Length") == "2"
@test result.resp.body == Vector{UInt8}("ok")

Routing.empty!()

end # module test_bukdu_head


text:
using MUMPSjInv
using Test
using LinearAlgebra
using SparseArrays
using Printf

include("getDivGrad.jl");
Ns = (8,16,24,32,48,64)
MUMPStime=zeros(length(Ns))
MUMPStimeOOC=zeros(length(Ns))
Juliatime=zeros(length(Ns))

for i=1:length(Ns)
	println(@sprintf("Solve div-grad on %d x %d x %d grid",Ns[i],Ns[i],Ns[i]));

	A = getDivGrad(Ns[i],Ns[i],Ns[i]);
	n = size(A,1);
    rhs = randn(n);

	# solve using MUMPSjInv
	MUMPStimeOOC[i] = @elapsed begin
		x = solveMUMPS(A,rhs,1,1);
	end
	# solve using MUMPSjInv in core
	MUMPStime[i] = @elapsed begin
		x = solveMUMPS(A,rhs,1,0);
	end

	# solve using MUMPSjInv
	Juliatime[i] = @elapsed begin
		x = A\rhs;
	end
end

for i=1:length(Ns)
	println(@sprintf("| %d^3 | %1.3f | %1.3f | %1.3f | %1.1f |",Ns[i],MUMPStimeOOC[i],MUMPStime[i],Juliatime[i],Juliatime[i]/MUMPStime[i]))
end


text:
# ---
# title: 667. Beautiful Arrangement II
# id: problem667
# author: Indigo
# date: 2021-06-28
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/beautiful-arrangement-ii/description/>
# hidden: true
# ---
#
# Given two integers `n` and `k`, you need to construct a list which contains
# `n` different positive integers ranging from `1` to `n` and obeys the
# following requirement:
# Suppose this list is [a1, a2, a3, ... , an], then the list [|a1 \- a2|, |a2 \-
# a3|, |a3 \- a4|, ... , |an-1 \- an|] has exactly `k` distinct integers.
#
# If there are multiple answers, print any of them.
#
# **Example 1:**
#
#
#
#     Input: n = 3, k = 1
#     Output: [1, 2, 3]
#     Explanation: The [1, 2, 3] has three different positive integers ranging from 1 to 3, and the [1, 1] has exactly 1 distinct integer: 1.
#
#
# **Example 2:**
#
#
#
#     Input: n = 3, k = 2
#     Output: [1, 3, 2]
#     Explanation: The [1, 3, 2] has three different positive integers ranging from 1 to 3, and the [2, 1] has exactly 2 distinct integers: 1 and 2.
#
#
# **Note:**
#
#   1. The `n` and `k` are in the range 1 <= k < n <= 104.
#
#
## @lc code=start
using LeetCode

function beautiful_arrangement(n, k)
    res = [1:(n-k);]
    i = res[end] + 1
    j = n
    while i < j
        push!(res, j)
        j -= 1
        push!(res, i)
        i += 1
    end
    i == j && push!(res, i)
    return res
end
## @lc code=end


text:
macro brenttrace()
    quote
        if tracing
            dt = Dict()
            if extended_trace
                dt["x_minimum"] = x_minimum
                dt["x_lower"] = x_lower
                dt["x_upper"] = x_upper
            end
            update!(tr,
                    iteration,
                    f_minimum,
                    NaN,
                    dt,
                    store_trace,
                    show_trace,
                    show_every,
                    callback)
        end
    end
end

immutable Brent <: Optimizer end

function optimize{T <: AbstractFloat}(
        f::Function, x_lower::T, x_upper::T,
        mo::Brent;
        rel_tol::T = sqrt(eps(T)),
        abs_tol::T = eps(T),
        iterations::Integer = 1_000,
        store_trace::Bool = false,
        show_trace::Bool = false,
        callback = nothing,
        show_every = 1,
        extended_trace::Bool = false)

    if !(x_lower < x_upper)
        error("x_lower must be less than x_upper")
    end

    # Save for later
    initial_lower = x_lower
    initial_upper = x_upper

    const golden_ratio::T = 0.5 * (3.0 - sqrt(5.0))

    x_minimum = x_lower + golden_ratio*(x_upper-x_lower)
    f_minimum = f(x_minimum)
    f_calls = 1 # Number of calls to f

    step = zero(T)
    step_old = zero(T)

    x_minimum_old = x_minimum
    x_minimum_old_old = x_minimum

    f_minimum_old = f_minimum
    f_minimum_old_old = f_minimum

    iteration = 0
    converged = false

    # Trace the history of states visited
    tr = OptimizationTrace{typeof(mo)}()
    tracing = store_trace || show_trace || extended_trace || callback != nothing
    @brenttrace

    while iteration < iterations

        p = zero(T)
        q = zero(T)

        tolx = rel_tol * abs(x_minimum) + abs_tol

        x_midpoint = (x_upper+x_lower)/2

        if abs(x_minimum - x_midpoint) <= 2*tolx - (x_upper-x_lower)/2
            converged = true
            break
        end

        iteration += 1

        if abs(step_old) > tolx
            # Compute parabola interpolation
            # x_minimum + p/q is the optimum of the parabola
            # Also, q is guaranteed to be positive

            r = (x_minimum - x_minimum_old) * (f_minimum - f_minimum_old_old)
            q = (x_minimum - x_minimum_old_old) * (f_minimum - f_minimum_old)
            p = (x_minimum - x_minimum_old_old) * q - (x_minimum - x_minimum_old) * r
            q = 2(q - r)

            if q > 0
                p = -p
            else
                q = -q
            end
        end

        if abs(p) < abs(q*step_old/2) && p < q*(x_upper-x_minimum) && p < q*(x_minimum-x_lower)
            step_old = step
            step = p/q

            # The function must not be evaluated too close to x_upper or x_lower
            x_temp = x_minimum + step
            if ((x_temp - x_lower) < 2*tolx || (x_upper - x_temp) < 2*tolx)
                step = (x_minimum < x_midpoint) ? tolx : -tolx
            end
        else
            step_old = (x_minimum < x_midpoint) ? x_upper - x_minimum : x_lower - x_minimum
            step = golden_ratio * step_old
        end

        # The function must not be evaluated too close to x_minimum
        if abs(step) >= tolx
            x_new = x_minimum + step
        else
            x_new = x_minimum + ((step > 0) ? tolx : -tolx)
        end

        f_new = f(x_new)
        f_calls += 1

        if f_new <= f_minimum
            if x_new < x_minimum
                x_upper = x_minimum
            else
                x_lower = x_minimum
            end
            x_minimum_old_old = x_minimum_old
            f_minimum_old_old = f_minimum_old
            x_minimum_old = x_minimum
            f_minimum_old = f_minimum
            x_minimum = x_new
            f_minimum = f_new
        else
            if x_new < x_minimum
                x_lower = x_new
            else
                x_upper = x_new
            end
            if f_new <= f_minimum_old || x_minimum_old == x_minimum
                x_minimum_old_old = x_minimum_old
                f_minimum_old_old = f_minimum_old
                x_minimum_old = x_new
                f_minimum_old = f_new
            elseif f_new <= f_minimum_old_old || x_minimum_old_old == x_minimum || x_minimum_old_old == x_minimum_old
                x_minimum_old_old = x_new
                f_minimum_old_old = f_new
            end
        end

        @brenttrace
    end

    return UnivariateOptimizationResults("Brent's Method",
                                         initial_lower,
                                         initial_upper,
                                         x_minimum,
                                         Float64(f_minimum),
                                         iteration,
                                         iteration == iterations,
                                         converged,
                                         rel_tol,
                                         abs_tol,
                                         tr,
                                         f_calls)
end


text:
function lsr(eigval_sol)
    eigval_sol = real.(eigval_sol)
    sort!(eigval_sol)
    gap = Array{Float64}(undef, length(eigval_sol) - 1)
    r = Array{Float64}(undef, length(eigval_sol) - 2)
    for i in 1:length(eigval_sol) - 1
        gap[i] = eigval_sol[i + 1] - eigval_sol[i]
    end

    for i in 1:length(eigval_sol) - 2
        r[i] = min(gap[i + 1], gap[i])/max(gap[i + 1], gap[i])
    end
    # gap_of_E
    mean(r)
end

function EE(eigvec_sol, left_chain, right_chain, nhil::Int)

    # entroy = Array{Float64}(undef, nhil)
    # for i in 1:nhil
    #     rho = sparse(left_chain, right_chain, eigvec_sol[:, i])
    #     rho_left = rho * rho'
    #     entroy[i] = log(tr(rho_left^2)) / (1 - 2)
    # end
    # mean(entroy[ceil(Int,nhil * 3/8):ceil(Int, nhil * 5/8)])

    # entroy = Array{Float64}(undef, nhil)
    # for i in 1:nhil
    #     rho = sparse(left_chain, right_chain, eigvec_sol[:, i])
    #     rho_left = rho * rho'
    #     rho_left = Array(rho_left)
    #     entroy[i] = real(-tr(rho_left * log(rho_left)))
    # end
    # mean(entroy)

    # entroy = Array{Float64}(undef, nhil)
    # # test = []
    # # my_flag::Bool = true
    # for i in 1:nhil
    #     rho = sparse(left_chain, right_chain, eigvec_sol[:, i])
    #     rho_left = rho * rho'
    #     rho_left = Symmetric(Array(rho_left))
    #     # eig_v = abs.(real(eigvals(rho_left)))
    #     eig_v = real.(eigvals(rho_left))
    #     # eig_v = LinearAlgebra.LAPACK.syev!('N', 'U', rho_left)
    #     entroy[i] = -sum(eig_v .* log.(abs.(eig_v)))
    #     if isnan(entroy[i])
    #         # push!(test, eig_v)
    #         entroy[i] = 0.0
    #     end
    #     # Symmetric[rho_left]
    #     # my_flag = my_flag && (sum(rho_left .== Symmetric(rho_left)) == length(rho_left))
    # end
    # mean(entroy[ceil(Int,nhil * 3/8):ceil(Int, nhil * 5/8)])    # 根据文献只取中间1/4
    # # test
    # # return my_flag

    entroy = Array{Float64}(undef, nhil)
    # test = []
    # my_flag::Bool = true
    for i in 1:nhil
        rho = sparse(left_chain, right_chain, eigvec_sol[:, i])

        #################################
        # GPU
        # rho = CuArray(Array(rho))
        # rho_left = rho * rho'
        # rho = nothing
        # # eig_v = CUSOLVER.syevd!( 'N', 'U', CuArray(Symmetric(Array(rho_left))) )
        # eig_v = CUSOLVER.syevd!( 'N', 'U', rho_left )
        # eig_v = Array(eig_v)
        #################################

        #################################
        # CPU
        rho_left = rho * rho'
        rho = nothing
        eig_v = eigvals( Symmetric(Array(rho_left)) )
        # eig_v = abs.(real(eigvals(rho_left)))
        # eig_v = real.(eigvals(rho_left))
        # eig_v = LinearAlgebra.LAPACK.syev!('N', 'U', rho_left)
        #################################

        entroy[i] = -sum(eig_v .* log.(abs.(eig_v)))
        if isnan(entroy[i])
            # push!(test, eig_v)
            entroy[i] = 0.0
        end
        # Symmetric[rho_left]
        # my_flag = my_flag && (sum(rho_left .== Symmetric(rho_left)) == length(rho_left))
    end
    mean(entroy[ceil(Int,nhil * 3/8):ceil(Int, nhil * 5/8)])    # 根据文献只取中间1/4
    # test
    # return my_flag
end