text:
"""
    Voltage(;name)

Acts as an ideal voltage source with no internal resistance.

# States:
See [OnePort](@ref)

# Connectors:
- `p` Positive pin
- `n` Negative pin
- `V` [RealInput](@ref) Input for the voltage control signal, i.e. `V ~ p.v - n.v`
"""
function Voltage(;name)
    @named oneport = OnePort()
    @unpack v, i = oneport
    @named V = RealInput()
    eqs = [
        v ~ V.u
    ]

    extend(ODESystem(eqs, t, [], []; name=name, systems=[V]), oneport)
end

"""
    Current(;name)

Acts as an ideal current source with no internal resistance.

# States:
See [OnePort](@ref)

# Connectors:
- `p` Positive pin
- `n` Negative pin
- `I` [RealInput](@ref) Input for the current control signal, i.e. `I ~ p.i
"""
function Current(;name)
    @named oneport = OnePort()
    @unpack v, i = oneport
    @named I = RealInput()
    eqs = [
        i ~ I.u
    ]

    extend(ODESystem(eqs, t, [], []; name=name, systems=[I]), oneport)
end


text:
cluster = [1 1 1 1; 2 2 3 3; 4 4 5 5]
X = Array(reshape(1:24, 2, 3, 4))

@testset "pool" begin
    @testset "GlobalPool" begin
        glb_cltr = [1 1 1 1; 1 1 1 1; 1 1 1 1]
        p = GlobalPool(:add, 3, 4)
        @test p(X) == sumpool(glb_cltr, X)
    end

    @testset "LocalPool" begin
        p = LocalPool(:add, cluster)
        @test p(X) == sumpool(cluster, X)
    end

    @testset "TopKPool" begin
        N = 10
        k, in_channel = 4, 7
        X = rand(in_channel, N)
        for T = [Bool, Float64]
            adj = rand(T, N, N)
            p = TopKPool(adj, k, in_channel)
            @test eltype(p.p) === Float32
            @test size(p.p) == (in_channel,)
            @test eltype(p.Ã) === T
            @test size(p.Ã) == (k, k)
            y = p(X)
            @test size(y) == (in_channel, k)
        end
    end
end

ext: jl


text:
# idx in A_aux matrix
# assuming that x, y can be in [1, 2, ... m*n+3]
function aux_matrix_connectivity_entries_to_lin_idx(x,y, m,n)
    if (x, y) == (1, 2)
      return 1
    elseif x == 2
      # y - 2 + 1
      return y - 1
    elseif y == m*n + 3
      #          y                - 2 + 1 ... for special case of right edge
      # (m*n + 3) + mod(x - 3, m) - 2 + 1
      return m*n + 2 + mod(x - 3, m)
    # general case within the matrix
    # ... horizontal ... indexed first
    elseif mod(x - y, m) == 0
      return y - 1
    # ... vertical
    else
      # max index for horiz. case  +  idx of x  -  col_number (because there are 1 less vertical connections than #of rows)
      # (m*n + 3) + m - 2          +  x - 2     -  div(x-3, n)
      return m*(n + 1) + x - 1 - div(x-3, m)
    end
end

#
function aux_matrix_connectivity_entries_to_lin_idx(x,y, m,n,s)
    if (x, y) == (1, 2)
      return 1
    elseif x == 2
      # left end of electrode -> horizontal
      # idx of y          + transpose to layer
      # y - 2 + 1         + m*(n+1)*div(y-3, m*n)
      return y - 1 + m*div(y-3, m*n)
    elseif y == m*n*s + 3
      # right end of electrode -> horizontal
      # 1 + #of standard horizontal of the first layer    + #row                + transpose to layer
      # 1 + m*n                                           + mod(x-3, m) + 1     + m*(n + 1)*div(x-3, m*n)
      return m*n + 2 + mod(x-3,m) + m*(n+1)*div(x-3, m*n)
    elseif div(x-3, m*n) == div(y-3, m*n)
      if mod(x - y, m) == 0
        # horizontal
        # idx of y          + transpose to layer
        # y - 2 + 1         + m*(n+1)*div(x-3, m*n)
        return y - 1 + m*div(y-3, m*n)
      else
        # vertical
        # 1 + #horizontal    + idx of x    - #col   (because there is one connection less than row entries in each column)
        # 1 + m*(n+1)*s     + x - 2       - div(x-3, m)
        return m*(n+1)*s - 1 + x - div(x-3, m)
      end
    else
      # interlayer
      # 1 + # horizontal  + #vertical     + idx of x
      # 1 + m*(n+1)*s     + (m-1)*n*s     + x - 2
      return (m*(2*n+1) - n)*s + x - 1
    end
end

function ans_322(in_f, x, y)
  f(x, y) = in_f(x, y, 3,2,2)

  if (x,y) == (1,2); return f(x,y) == 1
  elseif (x,y) == (2,3); return f(x,y) == 2
  elseif (x,y) == (2,4); return f(x,y) == 3
  elseif (x,y) == (2,5); return f(x,y) == 4
  #
  elseif (x,y) == (3,6); return f(x,y) == 5
  elseif (x,y) == (4,7); return f(x,y) == 6
  elseif (x,y) == (5,8); return f(x,y) == 7
  #
  elseif (x,y) == (6,15); return f(x,y) == 8
  elseif (x,y) == (7,15); return f(x,y) == 9
  elseif (x,y) == (8,15); return f(x,y) == 10
  #
  #
  elseif (x,y) == (2,9); return f(x,y) == 11
  elseif (x,y) == (2,10); return f(x,y) == 12
  elseif (x,y) == (2,11); return f(x,y) == 13
  #
  elseif (x,y) == (9,12); return f(x,y) == 14
  elseif (x,y) == (10,13); return f(x,y) == 15
  elseif (x,y) == (11,14); return f(x,y) == 16
  #
  elseif (x,y) == (12,15); return f(x,y) == 17
  elseif (x,y) == (13,15); return f(x,y) == 18
  elseif (x,y) == (14,15); return f(x,y) == 19
  else
    return "ach jo"
  end
end

max_lin_idx(m,n) =  if m==1
                      aux_matrix_connectivity_entries_to_lin_idx(m*n + 2, m*n + 3, m, n)
                    else
                      aux_matrix_connectivity_entries_to_lin_idx(m*n + 1, m*n + 2, m, n)
                    end
max_lin_idx(m,n,s) =  if s==1
                        max_lin_idx(m,n)
                      else
                        aux_matrix_connectivity_entries_to_lin_idx(m*n*(s-1) + 2, m*n*s + 2, m, n, s)
                      end

function add_value_to_matrix_row(matrix_header, equation_row, previous_id, next_id, value, dims)
    lin_idx = aux_matrix_connectivity_entries_to_lin_idx(previous_id, next_id, dims...)

    push!(matrix_header, (lin_idx, (previous_id, next_id)))

    push!(equation_row, (lin_idx, value))
end

get_previous_list(m, n) = [(0, -1), (-1, 0)]
get_next_list(m, n) = [(1, 0), (0, 1)]

get_previous_list(m, n, s) = [(0, -1, 0), (-1, 0, 0), (0, 0, -1)]
get_next_list(m, n, s) = [(1, 0, 0), (0, 1, 0), (0, 0, 1)]

function previous(i, j, k, auxilary_A)
  if j == 0
    return [1]
  else
    res = Array{Int64}(undef, length(get_previous_list(0,0,0)))
    for (l, d) in enumerate(get_previous_list(0,0,0))
      res[l] = auxilary_A[(i+1, j+1, k+1) .+ d...]
    end
    deleteat!(res, findall(x->x==-1,res))
    return res
  end
end

function previous(i, j, auxilary_A)
  if j == 0
    return [1]
  else
    res = Array{Int64}(undef, length(get_previous_list(0,0)))
    for (k, d) in enumerate(get_previous_list(0,0))
      res[k] = auxilary_A[(i+1, j+1) .+ d...]
    end
    deleteat!(res, findall(x->x==-1,res))
    return res
  end
end

function next(i,j,k, auxilary_A)
  if j == 0
    return auxilary_A[2:end-1, 2, 2:end-1][:]
  else
    res = Array{Int64}(undef, length(get_next_list(0,0,0)))
    for (l, d) in enumerate(get_next_list(0,0,0))
      res[l] = auxilary_A[(i+1, j+1, k+1) .+ d...]
    end
    deleteat!(res, findall(x->x==-1,res))
    return res
  end
end

function next(i,j, auxilary_A)
  if j == 0
    return auxilary_A[2:end-1, 2]
  else
    res = Array{Int64}(undef, length(get_next_list(0,0)))
    for (k, d) in enumerate(get_next_list(0,0))
      res[k] = auxilary_A[(i+1, j+1) .+ d...]
    end
    deleteat!(res, findall(x->x==-1,res))
    return res
  end
end

function add_row_to_sparse_input!(sparse_input,
                                  new_row,
                                  row_idx)
  for item in new_row
    push!(sparse_input[1], row_idx)
    push!(sparse_input[2], item[1])
    push!(sparse_input[3], typeof(item[2]) <: Number ? float(item[2]) : item[2])
  end
end

function add_equation_I_row(pos, auxilary_A, dims, matrix_header, sys_row_idx, sparse_input, RHS)
  new_row = []
  act_id = auxilary_A[pos .+ 1 ...]
  #
  for node_id in previous(pos..., auxilary_A)
      add_value_to_matrix_row(matrix_header, new_row, node_id, act_id, -1, dims)
  end
  for node_id in next(pos..., auxilary_A)
      add_value_to_matrix_row(matrix_header, new_row, act_id, node_id, 1, dims)
  end

  sort!(new_row, by = first)

  sys_row_idx[1] += 1
  add_row_to_sparse_input!(sparse_input, new_row, sys_row_idx[1])
  push!(RHS, 0.0)
  return
end

function add_U_eq_for_each_relevant_path!(Z_vector, auxilary_A::Array{<:Integer, 3}, dims, matrix_header, sys_row_idx, sparse_input, RHS)
  A = auxilary_A[2:end-1, 2:end-1, 2:end-1]
  m, n, s = dims

  for layer in 1:s
    for row in 1:m
      for fall_point in n+1:-1: (row != m ? 1 : n+1)
        new_row = []
        # first index of each path is 2
        idx1 = 2
        column = 1
        act_row = row
        while column < n + 1
          idx2 = A[act_row, column, layer]
          add_value_to_matrix_row(
              matrix_header,
              new_row,
              idx1, idx2,
              Z_vector[aux_matrix_connectivity_entries_to_lin_idx(idx1, idx2, dims...)],
              dims
          )
          idx1 = idx2
          #
          if fall_point == column && act_row == row
            act_row += 1
          else
            column += 1
          end
        end
        # last index of each path is m*n*s + 3
        idx2 = m*n*s + 3
        add_value_to_matrix_row(
            matrix_header,
            new_row,
            idx1, idx2,
            Z_vector[aux_matrix_connectivity_entries_to_lin_idx(idx1, idx2, dims...)],
            dims
        )
        #
        sort!(new_row, by = first)

        sys_row_idx[1] += 1
        add_row_to_sparse_input!(sparse_input, new_row, sys_row_idx[1])
        push!(RHS, 1.0)
      end
    end
  end

  for layer in 1:s-1
    for row in 1:m
      for fall_point in n:-1: 1
        new_row = []
        # first index of each path is 2
        idx1 = 2
        column = 1
        act_layer = layer
        while column < n + 1
          idx2 = A[row, column, act_layer]
          add_value_to_matrix_row(
              matrix_header,
              new_row,
              idx1, idx2,
              Z_vector[aux_matrix_connectivity_entries_to_lin_idx(idx1, idx2, dims...)],
              dims
          )
          idx1 = idx2
          #
          if fall_point == column && act_layer == layer
            act_layer += 1
          else
            column += 1
          end
        end
        # last index of each path is m*n*s + 3
        idx2 = m*n*s + 3
        add_value_to_matrix_row(
            matrix_header,
            new_row,
            idx1, idx2,
            Z_vector[aux_matrix_connectivity_entries_to_lin_idx(idx1, idx2, dims...)],
            dims
        )
        #
        sort!(new_row, by = first)

        sys_row_idx[1] += 1
        add_row_to_sparse_input!(sparse_input, new_row, sys_row_idx[1])
        push!(RHS, 1.0)
      end
    end
  end
  return
end

function add_U_eq_for_each_relevant_path!(Z_vector, auxilary_A::Array{<:Integer, 2}, dims, matrix_header, sys_row_idx, sparse_input, RHS)
  A = auxilary_A[2:end-1, 2:end-1]
  m,n = dims

  for row in 1:m
    for fall_point in n+1:-1: (row != m ? 1 : n+1)
      new_row = []
      # first index of each path is 2
      idx1 = 2
      column = 1
      act_row = row
      while column < n + 1
        idx2 = A[act_row, column]
        add_value_to_matrix_row(
            matrix_header,
            new_row,
            idx1, idx2,
            Z_vector[aux_matrix_connectivity_entries_to_lin_idx(idx1, idx2, dims...)],
            dims
        )
        idx1 = idx2
        #
        if fall_point == column && act_row == row
          act_row += 1
        else
          column += 1
        end
      end
      # last index of each path is m*n + 3
      idx2 = m*n + 3
      add_value_to_matrix_row(
          matrix_header,
          new_row,
          idx1, idx2,
          Z_vector[aux_matrix_connectivity_entries_to_lin_idx(idx1, idx2, dims...)],
          dims
      )
      #
      sort!(new_row, by = first)

      sys_row_idx[1] += 1
      add_row_to_sparse_input!(sparse_input, new_row, sys_row_idx[1])
      push!(RHS, 1.0)
    end
  end

  return
end

function get_readable_matrix_header(matrix_header)
  sort!(matrix_header, by = first)
  res = [matrix_header[1]]
  for item in matrix_header
    if item[1] != res[end][1]
      push!(res, item)
    end
  end
  return [item[2] for item in res]
end

# auxilary_A has 2 on left, m*n+3 on the right, -1 on top and bottom
function vector_to_lin_sys(Z_vector, auxilary_A)
  # sparse_input =  (     sparse_input_row_idxs,
  #                       sparse_input_col_idxs,
  #                       sparse_input_vals )
  matrix_header = []

  sparse_input = (Int64[], Int64[], Union{Float64, Function}[])
  sys_row_idx = [0]

  RHS = Float64[]

  dims = size(auxilary_A) .- 2

  mode_3D = length(dims) == 3
  # I verticesnext(i,j)
  # initial vertex "2"
  if mode_3D
    add_equation_I_row((1,0,1), auxilary_A, dims, matrix_header, sys_row_idx, sparse_input, RHS)
  else
    add_equation_I_row((1,0), auxilary_A, dims, matrix_header, sys_row_idx, sparse_input, RHS)
  end

  # other vertices
  if mode_3D
    for i in 1:dims[1], j in 1:dims[2], k in 1:dims[3]
      add_equation_I_row((i,j, k), auxilary_A, dims, matrix_header, sys_row_idx, sparse_input, RHS)
    end
  else
    for i in 1:dims[1], j in 1:dims[2]
        add_equation_I_row((i,j), auxilary_A, dims, matrix_header, sys_row_idx, sparse_input, RHS)
    end
  end

  # U equations
  add_U_eq_for_each_relevant_path!(Z_vector, auxilary_A, dims, matrix_header, sys_row_idx, sparse_input, RHS)

  matrix_header_output = get_readable_matrix_header(matrix_header)

  return matrix_header_output, sparse_input, RHS
end

function get_auxilary_graph_matrix(m, n)
    A_aux = Array{Int64}(undef, m+2, n+2)
    #
    aux_vec = Vector{Int64}(undef, m*n)
    for i in 1:length(aux_vec)
        aux_vec[i] = i+2
    end
    pre_A = reshape(aux_vec, m, n)
    #
    A_aux[2:end-1, 2:end-1] = pre_A
    A_aux[:, 1] .= 2
    A_aux[:, end] .= m*n+3
    A_aux[1, :] .= -1
    A_aux[end, :] .= -1
    return A_aux
end

function get_auxilary_graph_matrix(m, n, s)
    A_aux = Array{Int64}(undef, m+2, n+2, s+2)
    #
    aux_vec = Vector{Int64}(undef, m*n*s)
    for i in 1:length(aux_vec)
        aux_vec[i] = i+2
    end
    pre_A = reshape(aux_vec, m, n, s)
    #
    A_aux[2:end-1, 2:end-1, 2:end-1] = pre_A
    A_aux[:, 1, :] .= 2
    A_aux[:, end, :] .= m*n*s+3
    A_aux[1, :,  :] .= -1
    A_aux[end, :, :] .= -1
    A_aux[:, :, 1] .= -1
    A_aux[:, :, end] .= -1
    return A_aux
end

function get_large_material_matrix(pre_A::Array{<:Integer, 2})
  (m,n) = size(pre_A)
  A = Array{Int64}(undef, m+2, n+2)
  A[2:end-1, 2:end-1] = pre_A
  A[:, 1] .= i_LSM
  A[:, end] .= i_LSM
  A[1, :] .= -1
  A[end, :] .= -1
  return A
end

function get_large_material_matrix(pre_A::Array{<:Integer, 3})
    (m,n,s) = size(pre_A)
    A = Array{Int64}(undef, m+2, n+2, s+2)
    A[2:end-1, 2:end-1, 2:end-1] = pre_A
    A[:, 1, :] .= i_LSM
    A[:, end, :] .= i_LSM
    A[1, :, :] .= -1
    A[end, :, :] .= -1
    A[:, :, 1] .= -1
    A[:, :, end] .= -1
    return A
end

function add_Z_to_vector!(Z_vector, pos, dir, aux_A, large_material_A, p, dims)
    n1 = large_material_A[pos...]
    n2 = large_material_A[pos .+ dir...]
    if n2 == -1
        return
    end

    idx1 = aux_A[pos...]
    idx2 = aux_A[pos .+ dir...]
    if idx2 < idx1
        aux = idx1
        idx1 = idx2
        idx2 = aux
    end

    current_key = aux_matrix_connectivity_entries_to_lin_idx(idx1, idx2, dims...)
    if Z_vector[current_key] == 0
        Z_vector[current_key] =
            get_Z_entry_from_material_matrix_codes(n1, n2, p)
    else
        f = deepcopy(Z_vector[current_key])
        Z_vector[current_key] = w -> f(w) +
            get_Z_entry_from_material_matrix_codes(n1, n2, p)(w)
    end
end

function material_matrix_to_lin_sys(
            material_A::Array{<:Integer}=[0 0 0; 1 1 1],
            p::parameters=parameters()
    )

    dims = size(material_A)
    mode_3D = length(dims) == 3
    aux_A = get_auxilary_graph_matrix(dims...)
    large_material_A = get_large_material_matrix(material_A)

    Z_vector = Vector{Union{Function}}(undef, max_lin_idx(dims...))
    Z_vector .= w -> 0
    for d in vcat(get_previous_list(dims...), get_next_list(dims...))
      for i in 2:dims[1]+1
        for j in 2:dims[2]+1
          if mode_3D
            for k in 2:dims[3]+1
              add_Z_to_vector!(Z_vector, (i, j, k), d, aux_A, large_material_A, p, dims)
            end
          else
            add_Z_to_vector!(Z_vector, (i, j), d, aux_A, large_material_A, p, dims)
          end
        end
      end
    end
    for i in 2:dims[1]+1
      if mode_3D
        for k in 2:dims[3]+1
          add_Z_to_vector!(Z_vector, (i,1, k), (0, +1, 0), aux_A, large_material_A, p, dims)
          add_Z_to_vector!(Z_vector, (i,dims[2]+2, k), (0, -1, 0), aux_A, large_material_A, p, dims)
        end
      else
        add_Z_to_vector!(Z_vector, (i,1), (0, +1), aux_A, large_material_A, p, dims)
        add_Z_to_vector!(Z_vector, (i,dims[2]+2), (0, -1), aux_A, large_material_A, p, dims)
      end
    end
    #return Z_vector, aux_A, large_material_A
    return vector_to_lin_sys(Z_vector, aux_A)
end

max_stars_repo_name: Masicko/ImageToEIS

max_issues_repo_name: Masicko/ImageToEIS

max_forks_repo_name: Masicko/ImageToEIS


text:
using GAFramework
using GAFramework.CoordinateGA
using LinearAlgebra
using FileIO
using Random

# Solves for function a + 2b + 3c + 4d + 5e = 42
# by minimizing abs(a + 2b + 3c + 4d + 5e - 42)
m = FunctionModel(x -> abs(x[1] + 2x[2] + 3x[3] + 4x[4] + 5x[5] - 42),
    [-100.,-100.,-100.,-100.,-100], # minimum corner
    [100.,100,100.,100.,100]) # maximum corner in rectangle
# Prints fitness every iteration
state = GAState(m, ngen=100, npop=3_000, elite_fraction=0.01,
    params=Dict(:mutation_rate=>0.1, :print_fitness_iter=>10))
println(ga!(state))

# Solves for function a + 2b + 3c + 4d + 5e = 42
# by minimizing abs(a + 2b + 3c + 4d + 5e - 42)
# without restricting the result to the [-100,100] rectangle
# but still initializing the population in that rectangle
m = FunctionModel(x -> abs(x[1] + 2x[2] + 3x[3] + 4x[4] + 5x[5] - 42),
    [-100.,-100.,-100.,-100.,-100], # minimum corner
    [100.,100,100.,100.,100], false) # maximum corner in rectangle
# Prints fitness every iteration
state = GAState(m, ngen=100, npop=3_000, elite_fraction=0.01,
    params=Dict(:mutation_rate=>0.1, :print_fitness_iter=>10))
println(ga!(state))

# Minimizes exp(|x|)
m = FunctionModel(x -> exp(sqrt(dot(x,x))),
    [-1.,-1.,-1.,-1.,-1], # minimum corner
    [1.,1.,1.,1.,1]) # maximum corner in rectangle
# Saves creature to file every 10 iterations
state = GAState(m, ngen=100, npop=3_000, elite_fraction=0.01,
    params=Dict(:mutation_rate=>0.1, :print_fitness_iter=>10,
        :save_creature_iter=>10, :file_name_prefix=>"minexp_3000"))
println(ga!(state))

# Change maximum number of generations and then continue
state.ngen = 200
ga!(state)

# Minimizes |tan(x + 1000)|
m = FunctionModel(x -> abs(tanh(x[1]+1000)), [-2000.], [2000.])
# Saves full state to file every 10 iterations
state = GAState(m, ngen=100, npop=3_000, elite_fraction=0.01,
    params=Dict(:mutation_rate=>0.1, :print_fitness_iter=>10,
        :save_state_iter=>25, :file_name_prefix=>"tan1000_3000"))
println(ga!(state))

if false
    # Reload state from file and restart
    state = loadgastate("tan1000_3000_state_25.jld")

    ga!(state)
end
# It's possible to also change the other fields in state and then continue on
# but that is usually only when you made a mistake with a filename or something


text:
export HopfFibration
export HSVtoRGB

struct ⭕
    b::Array{Complex}
    f::Array{Float64,2}
    m::Array{Float64,4}
    c::Array{Float64,4}
    s::Integer
    r::Float64
    q::ℍ
    p::Array{Float64}
end

λ(z₁, z₂) = begin
    d = 1 .- imag.(z₂)
    y = [real.(z₁) ./ d imag.(z₁) ./ d real.(z₂) ./ d]
    c = similar(y)
    sigmoid(x) = 2 / (1 + exp(-x)) .- 1
    for i in 1:size(y, 1)
        c[i, :] = LinearAlgebra.normalize(y[i, :]) .* tanh(LinearAlgebra.norm(y[i, :]))
    end
    c
end

HSVtoRGB(hsv) = begin
    H, S, V = hsv
    C = V * S
    X = C * (1 - Base.abs((H / 60) % 2 - 1))
    m = V - C
    if 0 ≤ H < 60
        R′, G′, B′ = C, X, 0
    elseif 60 ≤ H < 120
        R′, G′, B′ = X, C, 0
    elseif 120 ≤ H < 180
        R′, G′, B′ = 0, C, X
    elseif 180 ≤ H < 240
        R′, G′, B′ = 0, X, C
    elseif 240 ≤ H < 300
        R′, G′, B′ = X, 0, C
    elseif 300 ≤ H < 360
        R′, G′, B′ = C, 0, X
    else
        R′, G′, B′ = rand(3)
    end
    R, G, B = R′ + m, G′ + m, B′ + m
    [R; G; B]
end

function get_fibers(b::Array{Complex},
                    f::Array{Float64,2},
                    s::Integer,
                    r::Float64,
                    q::ℍ,
                    offset::Array{Float64})
    samples = size(b, 1)
    θ₁, θ₂ = f[:, 1], f[:, 2]
    Q = [1.0; 0.0; 0.0; 0.0]
    nᵢ = [0; 0; 1]
    s2 = Integer(s ÷ 3)
    ψ = range(0, stop = 2pi, length = s2)
    zero = fill(0, s2)
    m = Array{Float64,4}(undef, samples, s, s2, 3)
    c = similar(m)
    construct(b, ξ₂, q) = begin
        samples = size(ξ₂, 1)
        z₁ = exp.(im .* (ξ₂ .+ imag.(b) ./ 2)) .* cos.(atan.(real.(b)))
        z₂ = exp.(im .* (ξ₂ .- imag.(b) ./ 2)) .* sin.(atan.(real.(b)))
        rotatedz₁ = Array{Complex}(undef, samples)
        rotatedz₂ =  similar(rotatedz₁)
        for i in 1:samples
            rotated = ℍ([real.(z₁[i]); imag.(z₁[i]); real.(z₂[i]); imag.(z₂[i])]) * q
            rotatedz₁[i] = Complex(rotated.q[1], rotated.q[2])
            rotatedz₂[i] = Complex(rotated.q[3], rotated.q[4])
        end
        rotatedz₁, rotatedz₂
    end
    for i in 1:samples
        ξ₂ = range(θ₁[i], stop = θ₂[i], length = s)
        z₁, z₂ = construct(b[i], ξ₂, q)
        p = λ(z₁, z₂)
        p′ = circshift(p, 1)
        P = [real(z₁[1]); imag(z₁[1]); real(z₂[1]); imag(z₂[1])]
        hue = acos(LinearAlgebra.dot(P, Q)) / pi
        rgb = HSVtoRGB([hue * 360; 1.0; 1.0])
        c[i, :, :, :] = reshape(repeat(rgb', s * s2), s, s2, 3)
        for j in 1:s
            n = LinearAlgebra.normalize(p′[j, :] - p[j, :])
            u = LinearAlgebra.normalize(LinearAlgebra.cross(nᵢ, n))
            β = acos(LinearAlgebra.dot(nᵢ, n)) / 2
            h2 = ℍ([cos(β); sin(β) .* u])
            circle = [r .* cos.(ψ) r .* sin.(ψ) zero]
            m[i, j, :, :] = 💞(circle,
                               h2) + repeat((p[j, :])', s2, 1) + repeat(offset', s2, 1)
        end
    end
    m, c
end

⭕(b::Array{Complex}, f::Array{Float64,2}, s::Int64, r::Float64, q::ℍ, p::Array{Float64}) = begin
    m, c = get_fibers(b, f, s, r, q, p)
    ⭕(b, f, m, c, s, r, q, p)
end


text:
using StrainRecon
using Base.Test

@testset "BCD" begin
M  = [1 1 1; 1 1 0; 1 0 1; 1 0 0; 0 1 1; 0 1 0; 0 0 1;]
(m,n) = size(M)
R     = StrainRecon.sdiag(ones(m)*100);
atrue = [.6; 0.3; 0.1]
Data  = M*atrue

paramGauss = getStrainReconParam(Data,n,R=R,pwr=2.0)
paramLaplace = getStrainReconParam(Data,n,R=R,pwr=1.0)

for param=[paramGauss,paramLaplace]
    @testset "testBCD pwr=$(param.pwr)" begin
        maxTry = 10
        A0   = zeros(n,maxTry)
        Abcd = zeros(n,maxTry)
        Mbcd =  falses(m,n,maxTry)
        Jbcd = zeros(maxTry)
        mbcd = falses(m,n)
        for k=1:maxTry
            a0 = randomw(n)
            mbcd,abcd,his = bcd(param, a0,out=-1)
            Abcd[:,k] = abcd
        	Mbcd[:,:,k] = mbcd
            A0[:,k] = a0
            Jbcd[k] = misfit(param,mbcd,abcd)[1]
        end
        # find global minimizer
        kbcd = indmin(Jbcd)
        abcd = Abcd[:,kbcd]
        Mbcd = Mbcd[:,:,kbcd]

        @test norm(sum(Abcd,1)-1,Inf) < 1e-8
        @test norm(abcd - atrue)/norm(atrue) < 1e-8
        @test all(Abcd.>=0)
        @test all(Mbcd.==M)
        @test norm(Mbcd*abcd-Data)/norm(Data) < 1e-8
    end
end
end


text:
"""
NOTE: in julia we can represent closed abstract types with two different systems:
1. either Union{Type1, Type2}
2. or with a standard abstract type, where all type definitions are on the abstract type level
  (if you would only define functions on the concretetypes, the the inference engine will do weird things
   as it assumes there might be another concretetype in the future)

Approach 1.
===========

Advantages: Providing definitions for the concrete is actually enough for typeinference. More code reuse.
Disadvantages: The typeinference like for `Dict(:a => Some(1), :b => None())` would infer
  `Dict{Symbol, Any}` instead of `Dict{Symbol, Union{Some, None}}`.

  The internal details for this is that in many cases `Base.promote_typejoin(type1, type2)` is used to come up
  with a common type.
  While the type hierarchy is used for this, `Base.promote_typejoin` of two unrelated types will always result in `Any`.

  You would need to overload `promote_typejoin` to support
  `Dict(:a => Identity(1), :b => nothing) isa Dict{Symbol, Option{Int}}`.

Approach 2.
===========

Advantages: `Dict(:a => Some(1), :b => None())` would indeed infer `Dict{Symbol, Option}`
Disadvantages: you need to be careful to always implement functionalities first on separate functions unique to the
  sealed type, and then point generic functions to the specific one via `genericfunction(o::Option) = optionfunction(o)`

As we managed to make the promote_typejoin work as intended for `Union`, we follow Approach 1, in order to enable way
higher code-reuse and datatype-reuse.
"""
module DataTypesBasic

export Identity, isidentity, Const,
  Option, isoption, issome, isnone, iftrue, iffalse, getOption,
  Either, either, @either, iseither, isleft, isright, getleft, getright, getleftOption, getrightOption, flip_left_right,
  Try, @Try, @TryCatch, istry, issuccess, isfailure, Thrown, MultipleExceptions,
  ContextManager, @ContextManager

using Compat

# type definitions
include("Identity.jl")
include("Const.jl")
include("Option.jl")
include("Try.jl")
include("Either.jl")
include("ContextManager.jl")

# interoperability between types
include("promote_type.jl")
include("convert.jl")

end # module


text:
function testentropy()

  # check my implementation of entropy against that of Distributions.jl

  D = 5
  A = randn(D, D)
  Σ = A*A'
  μ = randn(D)

  @test Distributions.entropy(MvNormal(μ, Σ)) ≈ ApproximateVI.entropy(A)

  D = 11
  A = randn(D)
  Σ = Diagonal(A.^2)
  μ = randn(D)

  @test Distributions.entropy(MvNormal(μ, Σ)) ≈ ApproximateVI.entropy(A)

  σ = rand()+1e-4
  μ = randn()

  @test Distributions.entropy(Normal(μ, σ)) ≈ ApproximateVI.entropy(σ)

end


text:
using Base.Iterators: take, zip
using Test
using Arbitrary
using VectorSpaces

I2 = Vec{2, Int}
@test VectorSpaces.isvectorspace(I2)
F2 = Vec{2, Float64}
@test VectorSpaces.isvectorspace(F2)
F3 = Vec{3, Float64}
@test VectorSpaces.isvectorspace(F3)

F5 = VProd{F2, F3, Float64}
@test VectorSpaces.isvectorspace(F5)

F0 = VEmpty{Float64}
@test VectorSpaces.isvectorspace(F0)
F0′ = VProd{F0, F0, Float64}
@test VectorSpaces.isvectorspace(F0′)
F3′ = VProd{F3, F0, Float64}
@test VectorSpaces.isvectorspace(F3′)
F10 = VProd{F5, F2, Float64}
@test VectorSpaces.isvectorspace(F10)

F6 = VComp{Vec{2, Vec{3, Float64}}, Float64}
@test VectorSpaces.isvectorspace(F6)

F1 = VUnit{Float64}
@test VectorSpaces.isvectorspace(F1)
F1′ = VComp{VUnit{VUnit{Float64}}, Float64}
@test VectorSpaces.isvectorspace(F1′)
F3′′ = VComp{Vec{3, VUnit{Float64}}, Float64}
@test VectorSpaces.isvectorspace(F3′′)
F6′ = VComp{retype(F6, F1), Float64}
@test VectorSpaces.isvectorspace(F6′)

const BigRational = Rational{BigInt}

R = BigRational
CR = VComplex{R}
QR = VQuaternion{R}
OR = VOctonion{R}
SR = VSedenion{R}
@test VectorSpaces.isvectorspace(CR)
@test VectorSpaces.isvectorspace(QR)
@test VectorSpaces.isvectorspace(OR)
@test VectorSpaces.isvectorspace(SR)

@testset "Scalar" begin
    z = zero(R)
    e = one(R)

    as = arbitrary(R)
    bs = arbitrary(R)
    cs = arbitrary(R)

    for (a, b, c) in take(zip(as, bs, cs), 100)

        # Commutativity
        @test isequal(a + b, b + a)
        # Associativity
        @test isequal((a + b) + c, a + (b + c))
        # Identity
        @test isequal(z + a, a)
        @test isequal(a + z, a)
        # Inverse
        @test isequal((-a) + a, z)
        @test isequal(a + (-a), z)
        @test isequal(-(-a), a)
        # Subtraction
        @test isequal(a - b, a + (-b))

        # Commutativity
        @test isequal(a * b, b * a)
        # Associativity
        @test isequal((a * b) * c, a * (b * c))
        # Identity
        @test isequal(e * a, a)
        @test isequal(a * e, a)
        # Inverse
        if !isequal(a, z)
            @test isequal(inv(a) * a, e)
            @test isequal(a * inv(a), e)
            @test isequal(inv(inv(a)), a)
        end
        # Division
        if !isequal(b, z)
            @test isequal(a / b, a * inv(b))
        end

        # Distributivity
        @test isequal(a * (b + c), a * b + a * c)
        @test isequal((a + b) * c, a * c + b * c)

        # Conjugate
        @test isequal(conj(conj(a)), a)
        isreal(a) = isequal(a, conj(a))
        @test isreal(conj(a) * a)
    end
end

@testset "Complex" begin
    z = zero(CR)
    e = one(CR)

    as = arbitrary(CR)
    bs = arbitrary(CR)
    cs = arbitrary(CR)

    for (a, b, c) in take(zip(as, bs, cs), 100)

        @test length(a) == 2

        # Commutativity
        @test isequal(a + b, b + a)
        # Associativity
        @test isequal((a + b) + c, a + (b + c))
        # Identity
        @test isequal(z + a, a)
        @test isequal(a + z, a)
        # Inverse
        @test isequal((-a) + a, z)
        @test isequal(a + (-a), z)
        @test isequal(-(-a), a)
        # Subtraction
        @test isequal(a - b, a + (-b))

        # Commutativity
        @test isequal(a * b, b * a)
        # Associativity
        @test isequal((a * b) * c, a * (b * c))
        # Identity
        @test isequal(e * a, a)
        @test isequal(a * e, a)
        # Inverse
        if !isequal(a, z)
            @test isequal(inv(a) * a, e)
            @test isequal(a * inv(a), e)
            @test isequal(inv(inv(a)), a)
        end
        # Division
        if !isequal(b, z)
            @test isequal(a / b, a * inv(b))
        end

        # Distributivity
        @test isequal(a * (b + c), a * b + a * c)
        @test isequal((a + b) * c, a * c + b * c)

        # Conjugate
        @test isequal(conj(conj(a)), a)
        isreal(a) = isequal(a, conj(a))
        @test isreal(conj(a) * a)
    end
end

@testset "Quaternion" begin
    z = zero(QR)
    e = one(QR)

    as = arbitrary(QR)
    bs = arbitrary(QR)
    cs = arbitrary(QR)

    for (a, b, c) in take(zip(as, bs, cs), 100)

        @test length(a) == 4

        # Commutativity
        @test isequal(a + b, b + a)
        # Associativity
        @test isequal((a + b) + c, a + (b + c))
        # Identity
        @test isequal(z + a, a)
        @test isequal(a + z, a)
        # Inverse
        @test isequal((-a) + a, z)
        @test isequal(a + (-a), z)
        @test isequal(-(-a), a)
        # Subtraction
        @test isequal(a - b, a + (-b))

        # Associativity
        @test isequal((a * b) * c, a * (b * c))
        # Identity
        @test isequal(e * a, a)
        @test isequal(a * e, a)
        # Inverse
        if !isequal(a, z)
            @test isequal(inv(a) * a, e)
            @test isequal(a * inv(a), e)
            @test isequal(inv(inv(a)), a)
        end
        # Division
        if !isequal(b, z)
            @test isequal(a / b, a * inv(b))
        end

        # Distributivity
        @test isequal(a * (b + c), a * b + a * c)
        @test isequal((a + b) * c, a * c + b * c)

        # Conjugate
        @test isequal(conj(conj(a)), a)
        isreal(a) = isequal(a, conj(a))
        @test isreal(conj(a) * a)
    end
end

@testset "Octonion" begin
    z = zero(OR)
    e = one(OR)

    as = arbitrary(OR)
    bs = arbitrary(OR)
    cs = arbitrary(OR)

    for (a, b, c) in take(zip(as, bs, cs), 100)

        @test length(a) == 8

        # Commutativity
        @test isequal(a + b, b + a)
        # Associativity
        @test isequal((a + b) + c, a + (b + c))
        # Identity
        @test isequal(z + a, a)
        @test isequal(a + z, a)
        # Inverse
        @test isequal((-a) + a, z)
        @test isequal(a + (-a), z)
        @test isequal(-(-a), a)
        # Subtraction
        @test isequal(a - b, a + (-b))

        # Alternativity
        @test isequal((a * a) * b, a * (a * b))
        @test isequal((a * b) * b, a * (b * b))
        # Identity
        @test isequal(e * a, a)
        @test isequal(a * e, a)
        # Inverse
        if !isequal(a, z)
            @test isequal(inv(a) * a, e)
            @test isequal(a * inv(a), e)
            @test isequal(inv(inv(a)), a)
        end
        # Division
        if !isequal(b, z)
            @test isequal(a / b, a * inv(b))
        end

        # Distributivity
        @test isequal(a * (b + c), a * b + a * c)
        @test isequal((a + b) * c, a * c + b * c)

        # Conjugate
        @test isequal(conj(conj(a)), a)
        isreal(a) = isequal(a, conj(a))
        @test isreal(conj(a) * a)
    end
end

@testset "Sedenion" begin
    z = zero(SR)
    e = one(SR)

    as = arbitrary(SR)
    bs = arbitrary(SR)
    cs = arbitrary(SR)

    for (a, b, c) in take(zip(as, bs, cs), 100)

        @test length(a) == 16

        # Commutativity
        @test isequal(a + b, b + a)
        # Associativity
        @test isequal((a + b) + c, a + (b + c))
        # Identity
        @test isequal(z + a, a)
        @test isequal(a + z, a)
        # Inverse
        @test isequal((-a) + a, z)
        @test isequal(a + (-a), z)
        @test isequal(-(-a), a)
        # Subtraction
        @test isequal(a - b, a + (-b))

        # Power associativity
        @test isequal((a * a) * a, a * (a * a))
        # Identity
        @test isequal(e * a, a)
        @test isequal(a * e, a)
        # Inverse
        if !isequal(a, z)
            @test isequal(inv(a) * a, e)
            @test isequal(a * inv(a), e)
            @test isequal(inv(inv(a)), a)
        end
        # Division
        if !isequal(b, z)
            @test isequal(a / b, a * inv(b))
        end

        # Distributivity
        @test isequal(a * (b + c), a * b + a * c)
        @test isequal((a + b) * c, a * c + b * c)

        # Conjugate
        @test isequal(conj(conj(a)), a)
        isreal(a) = isequal(a, conj(a))
        @test isreal(conj(a) * a)
    end
end

Dmax = 4
for D in 0:Dmax, N in 0:D
    @testset "Forms D=$D N=$N" begin
        T = typeof(Form{D, N, BigInt}())
        z = zero(T)

        ntests = min(20, 1000 ÷ length(z))

        as = arbitrary(T)
        bs = arbitrary(T)
        cs = arbitrary(T)

        for (a, b, c) in take(zip(as, bs, cs), ntests)

            # Commutativity
            @test isequal(a + b, b + a)
            # Associativity
            @test isequal((a + b) + c, a + (b + c))
            # Identity
            @test isequal(z + a, a)
            @test isequal(a + z, a)
            # Inverse
            @test isequal((-a) + a, z)
            @test isequal(a + (-a), z)
            @test isequal(-(-a), a)
            # Subtraction
            @test isequal(a - b, a + (-b))

        end

        S = Rational{BigInt}
        e = Form{D, 0}(Vec((S(1),)))

        begin
            N1 = N
            T1 = typeof(Form{D, N1, S}())

            as = arbitrary(T1)

            for a in take(as, ntests)

                # Unit element
                @test isequal(e ∧ a, a)
                @test isequal(a ∧ e, a)

            end
        end

        for N1 in 0:N
            N2 = N - N1
            T1 = typeof(Form{D, N1, S}())
            T2 = typeof(Form{D, N2, S}())

            as = arbitrary(T1)
            bs = arbitrary(T2)
            ss = arbitrary(S)

            for (a, b, s) in take(zip(as, bs, ss), ntests)
                @test (a ∧ b) isa Form{D, N, S}

                # Linearity
                @test isequal((s * a) ∧ b, s * (a ∧ b))
                @test isequal(a ∧ (s * b), s * (a ∧ b))

                # Zero element
                @test isequal(zero(T1) ∧ b, zero(T))
                @test isequal(a ∧ zero(T2), zero(T))

                # (Anti-)Commutativity
                s = S((-1) ^ (N1 * N2))
                @test isequal(a ∧ b, s * (b ∧ a))

            end
        end

        for N1 in 0:N, N2 in 0:N-N1
            N3 = N - N1 - N2
            T1 = typeof(Form{D, N1, S}())
            T2 = typeof(Form{D, N2, S}())
            T3 = typeof(Form{D, N3, S}())

            as = arbitrary(T1)
            bs = arbitrary(T2)
            cs = arbitrary(T3)

            for (a, b, c) in take(zip(as, bs, cs), ntests)

                # Associativity
                @test isequal((a ∧ b) ∧ c, a ∧ (b ∧ c))

            end
        end

    end
end


text:
import LowRankModels: sort_observations
import FirstOrderOptimization: HopefulStepSize
using LowRankModels, FullRankModels, Compat
using Convex

srand(1)
# parameters
m, n, k = 100, 100, 3

# generate data
X0 = randn(m,k)
Y0 = randn(k,n)
A = X0*Y0 #+ .1*randn(m,n)
u,s,v = svd(A)
tau = sum(s)
P = sparse(float(randn(m,n) .>= .1)) # observed entries

I,J = findn(P) # observed indices (vectors)
@compat obs = Tuple{Int,Int}[(I[a],J[a]) for a = 1:length(I)]
observed_features, observed_examples = sort_observations(obs,size(P)...)

losses = fill(QuadLoss(1), n)
reg = TraceNormConstraint(tau)
gfrm = GFRM(A, losses, reg, k, observed_features, observed_examples,
	zeros(m,n), zeros(m+n,m+n))

params = FrankWolfeParams(maxiters = 100, reltol = 1e-10, stepsize = HopefulStepSize(1, 1, .5, 1.2, .01));
@time X_sketched, ch = fit!(gfrm, params)
@time X_sketched, ch = fit!(gfrm, params)

error()

Xhat = Array(X_sketched)

# compare with SDP solver
U = Variable(m,n)
W = Variable(m+n, m+n)

obj = sumsquares(P.*(U-A))
p = minimize(obj,
	nuclearnorm(U) <= tau)
solve!(p)

println("mean square error is ", vecnorm(U.value - Xhat) / vecnorm(U.value))
U.value = Xhat
println("objective of Convex problem evaluated at frank wolfe solution is $(Convex.evaluate(obj))
")


text:
using AutoGrad, Test
using AutoGrad: gradcheck, randcheck, gcheck, @gcheck


text:
"All converter variables"
function variable_dc_converter(pm::AbstractPowerModel; kwargs...)
    variable_conv_tranformer_flow(pm; kwargs...)
    variable_conv_reactor_flow(pm; kwargs...)

    variable_converter_active_power(pm; kwargs...)
    variable_converter_reactive_power(pm; kwargs...)
    variable_acside_current(pm; kwargs...)
    variable_dcside_power(pm; kwargs...)
    variable_converter_firing_angle(pm; kwargs...)

    variable_converter_filter_voltage(pm; kwargs...)
    variable_converter_internal_voltage(pm; kwargs...)

    variable_converter_to_grid_active_power(pm; kwargs...)
    variable_converter_to_grid_reactive_power(pm; kwargs...)
end

function variable_conv_tranformer_flow(pm::AbstractPowerModel; kwargs...)
    variable_conv_transformer_active_power_to(pm; kwargs...)
    variable_conv_transformer_reactive_power_to(pm; kwargs...)
end

function variable_conv_reactor_flow(pm::AbstractPowerModel; kwargs...)
    variable_conv_reactor_active_power_from(pm; kwargs...)
    variable_conv_reactor_reactive_power_from(pm; kwargs...)
end

"variable: `pconv_ac[j]` for `j` in `convdc`"
function variable_converter_active_power(pm::AbstractPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded::Bool = true)
    if bounded
        PowerModels.var(pm, nw, cnd)[:pconv_ac] = @variable(pm.model,
        [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_pconv_ac",
        lower_bound = PowerModels.ref(pm, nw, :convdc, i, "Pacmin", cnd),
        upper_bound = PowerModels.ref(pm, nw, :convdc, i, "Pacmax", cnd),
        start = PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "P_g", cnd, 1.0)
        )
    else
        PowerModels.var(pm, nw, cnd)[:pconv_ac] = @variable(pm.model,
        [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_pconv_ac",
        start = PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "P_g", cnd, 1.0)
        )
    end
end

"variable: `qconv_ac[j]` for `j` in `convdc`"
function variable_converter_reactive_power(pm::AbstractPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded::Bool = true)
    if bounded
        PowerModels.var(pm, nw, cnd)[:qconv_ac] = @variable(pm.model,
        [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_qconv_ac",
        lower_bound = PowerModels.ref(pm, nw, :convdc, i, "Qacmin", cnd),
        upper_bound = PowerModels.ref(pm, nw, :convdc, i, "Qacmax", cnd),
        start = PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "Q_g", cnd, 1.0)
        )
    else
        PowerModels.var(pm, nw, cnd)[:qconv_ac] = @variable(pm.model,
        [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_qconv_ac",
        start = PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "Q_g", cnd, 1.0)
        )
    end
end

"variable: `pconv_grid_ac_to[j]` for `j` in `convdc`"
function variable_conv_transformer_active_power_to(pm::AbstractPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded::Bool = true)
    bigM = 2;
    if bounded
        PowerModels.var(pm, nw, cnd)[:pconv_tf_to] = @variable(pm.model,
        [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_pconv_tf_to",
        lower_bound = -PowerModels.ref(pm, nw, :convdc, i, "Pacrated", cnd) * bigM,
        upper_bound =  PowerModels.ref(pm, nw, :convdc, i, "Pacrated", cnd) * bigM,
        start = PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "P_g", cnd, 1.0)
        )
    else
        PowerModels.var(pm, nw, cnd)[:pconv_tf_to] = @variable(pm.model,
        [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_pconv_tf_to",
        start = PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "P_g", cnd, 1.0)
        )
    end
end

"variable: `qconv_grid_ac_to[j]` for `j` in `convdc`"
function variable_conv_transformer_reactive_power_to(pm::AbstractPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded::Bool = true)
    bigM = 2;
    if bounded
        PowerModels.var(pm, nw, cnd)[:qconv_tf_to] = @variable(pm.model,
        [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_qconv_tf_to",
        lower_bound = -PowerModels.ref(pm, nw, :convdc, i, "Qacrated", cnd) * bigM,
        upper_bound =  PowerModels.ref(pm, nw, :convdc, i, "Qacrated", cnd) * bigM,
        start = PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "Q_g", cnd, 1.0)
        )
    else
        PowerModels.var(pm, nw, cnd)[:qconv_tf_to] = @variable(pm.model,
        [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_qconv_tf_to",
        start = PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "Q_g", cnd, 1.0)
        )
    end
end

"variable: `pconv_pr_from[j]` for `j` in `convdc`"
function variable_conv_reactor_active_power_from(pm::AbstractPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded::Bool = true)
    bigM = 2;
    if bounded
        PowerModels.var(pm, nw, cnd)[:pconv_pr_fr] = @variable(pm.model,
        [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_pconv_pr_fr",
        lower_bound = -PowerModels.ref(pm, nw, :convdc, i, "Pacrated", cnd) * bigM,
        upper_bound =  PowerModels.ref(pm, nw, :convdc, i, "Pacrated", cnd) * bigM,
        start = PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "P_g", cnd, 1.0)
        )
    else
        PowerModels.var(pm, nw, cnd)[:pconv_pr_fr] = @variable(pm.model,
        [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_pconv_pr_fr",
        start = PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "P_g", cnd, 1.0)
        )
    end
end

"variable: `qconv_pr_from[j]` for `j` in `convdc`"
function variable_conv_reactor_reactive_power_from(pm::AbstractPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded::Bool = true)
    bigM = 2;
    if bounded
        PowerModels.var(pm, nw, cnd)[:qconv_pr_fr] = @variable(pm.model,
        [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_qconv_pr_fr",
        lower_bound = -PowerModels.ref(pm, nw, :convdc, i, "Qacrated", cnd) * bigM,
        upper_bound =  PowerModels.ref(pm, nw, :convdc, i, "Qacrated", cnd) * bigM,
        start = PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "Q_g", cnd, 1.0)
        )
    else
        PowerModels.var(pm, nw, cnd)[:qconv_pr_fr] = @variable(pm.model,
        [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_qconv_pr_fr",
        start = PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "Q_g", cnd, 1.0)
        )
    end
end

"variable: `pconv_grid_ac[j]` for `j` in `convdc`"
function variable_converter_to_grid_active_power(pm::AbstractPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded::Bool = true)
    bigM = 2;
    if bounded
        PowerModels.var(pm, nw, cnd)[:pconv_tf_fr] = @variable(pm.model,
        [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_pconv_tf_fr",
        lower_bound = -PowerModels.ref(pm, nw, :convdc, i, "Pacrated", cnd) * bigM,
        upper_bound =  PowerModels.ref(pm, nw, :convdc, i, "Pacrated", cnd) * bigM,
        start = PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "P_g", cnd, 1.0)
        )
    else
        PowerModels.var(pm, nw, cnd)[:pconv_tf_fr] = @variable(pm.model,
        [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_pconv_tf_fr",
        start = PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "P_g", cnd, 1.0)
        )
    end
end

"variable: `qconv_grid_ac[j]` for `j` in `convdc`"
function variable_converter_to_grid_reactive_power(pm::AbstractPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded::Bool = true)
    bigM = 2;
    if bounded
        PowerModels.var(pm, nw, cnd)[:qconv_tf_fr] = @variable(pm.model,
        [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_qconv_tf_fr",
        lower_bound = -PowerModels.ref(pm, nw, :convdc, i, "Qacrated", cnd) * bigM,
        upper_bound =  PowerModels.ref(pm, nw, :convdc, i, "Qacrated", cnd) * bigM,
        start = PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "Q_g", cnd, 1.0)
        )
    else
        PowerModels.var(pm, nw, cnd)[:qconv_tf_fr] = @variable(pm.model,
        [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_qconv_tf_fr",
        start = PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "Q_g", cnd, 1.0)
        )
    end
end

"variable: `pconv_dc[j]` for `j` in `convdc`"
function variable_dcside_power(pm::AbstractPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded::Bool = true)
    bigM = 1.2; # to account for losses, maximum losses to be derived
    if bounded
        PowerModels.var(pm, nw, cnd)[:pconv_dc] = @variable(pm.model,
        [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_pconv_dc",
        lower_bound = -PowerModels.ref(pm, nw, :convdc, i, "Pacrated", cnd) * bigM,
        upper_bound =  PowerModels.ref(pm, nw, :convdc, i, "Pacrated", cnd) * bigM,
        start = PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "Pdcset", cnd, 1.0)
        )
    else
        PowerModels.var(pm, nw, cnd)[:pconv_dc] = @variable(pm.model,
        [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_pconv_dc",
        start = PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "Pdcset", cnd, 1.0)
        )
    end
end

"variable: `pconv_dc[j]` for `j` in `convdc`"
function variable_converter_firing_angle(pm::AbstractPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded::Bool = true)
    if bounded
        PowerModels.var(pm, nw, cnd)[:phiconv] = @variable(pm.model,
        [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_phiconv",
        lower_bound = 0,
        upper_bound = pi,
        start = acos(PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "Pdcset", cnd, 1.0) / sqrt((PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "Pacrated", cnd, 1.0))^2 + (PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "Qacrated", cnd, 1.0))^2))
        )
    else
        PowerModels.var(pm, nw, cnd)[:phiconv] = @variable(pm.model,
        [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_phiconv",
        start = acos(PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "Pdcset", cnd, 1.0) / sqrt((PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "Pacrated", cnd, 1.0))^2 + (PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "Qacrated", cnd, 1.0))^2))
        )
    end
end

"variable: `iconv_ac[j]` for `j` in `convdc`"
function variable_acside_current(pm::AbstractPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded::Bool = true)
    PowerModels.var(pm, nw, cnd)[:iconv_ac] = @variable(pm.model,
    [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_iconv_ac",
    lower_bound = 0,
    upper_bound = PowerModels.ref(pm, nw, :convdc, i, "Imax", cnd),
    start = PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "P_g", cnd, 1.0)
    )
end

"variable: `iconv_ac[j]` and `iconv_ac_sq[j]` for `j` in `convdc`"
function variable_acside_current(pm::PowerModels.AbstractWModels; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded::Bool = true)
    PowerModels.var(pm, nw, cnd)[:iconv_ac] = @variable(pm.model,
    [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_iconv_ac",
    lower_bound = 0,
    upper_bound = PowerModels.ref(pm, nw, :convdc, i, "Imax", cnd),
    start = PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "P_g", cnd, 1.0)
    )
    PowerModels.var(pm, nw, cnd)[:iconv_ac_sq] = @variable(pm.model,
    [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_iconv_ac_sq",
    lower_bound = 0,
    upper_bound = PowerModels.ref(pm, nw, :convdc, i, "Imax", cnd)^2,
    start = PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "P_g", cnd, 1.0)
    )
end

"variable: `itf_sq[j]` for `j` in `convdc`"
function variable_conv_transformer_current_sqr(pm::PowerModels.AbstractWModels; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded::Bool = true)
    bigM = 2; #TODO derive exact bound
    PowerModels.var(pm, nw, cnd)[:itf_sq] = @variable(pm.model,
    [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_itf_sq",
    lower_bound = 0,
    upper_bound = (bigM * PowerModels.ref(pm, nw, :convdc, i, "Imax", cnd))^2,
    start = PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "P_g", cnd, 1.0)^2
    )
end

"variable: `irc_sq[j]` for `j` in `convdc`"
function variable_conv_reactor_current_sqr(pm::PowerModels.AbstractWModels; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded::Bool = true)
    bigM = 2; #TODO derive exact bound
    PowerModels.var(pm, nw, cnd)[:irc_sq] = @variable(pm.model,
    [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_irc_sq",
    lower_bound = 0,
    upper_bound = (bigM * PowerModels.ref(pm, nw, :convdc, i, "Imax", cnd))^2,
    start = PowerModels.comp_start_value(ref(pm, nw, :convdc, i), "P_g", cnd, 1.0)^2
    )
end

function variable_converter_filter_voltage(pm::AbstractPowerModel; kwargs...)
    variable_converter_filter_voltage_magnitude(pm; kwargs...)
    variable_converter_filter_voltage_angle(pm; kwargs...)
end

"variable: `vmf[j]` for `j` in `convdc`"
function variable_converter_filter_voltage_magnitude(pm::AbstractPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded::Bool = true)
    bigM = 1.2; # only internal converter voltage is strictly regulated
    PowerModels.var(pm, nw, cnd)[:vmf] = @variable(pm.model,
    [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_vmf",
    lower_bound = PowerModels.ref(pm, nw, :convdc, i, "Vmmin", cnd) / bigM,
    upper_bound = PowerModels.ref(pm, nw, :convdc, i, "Vmmax", cnd) * bigM,
    start = PowerModels.ref(pm, nw, :convdc, i, "Vtar", cnd)
    )
end

"variable: `vaf[j]` for `j` in `convdc`"
function variable_converter_filter_voltage_angle(pm::AbstractPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded::Bool = true)
    bigM = 2*pi; #
    PowerModels.var(pm, nw, cnd)[:vaf] = @variable(pm.model,
    [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_vaf",
    lower_bound = -bigM,
    upper_bound =  bigM,
    start = 0
    )
end

function variable_converter_internal_voltage(pm::AbstractPowerModel; kwargs...)
    variable_converter_internal_voltage_magnitude(pm; kwargs...)
    variable_converter_internal_voltage_angle(pm; kwargs...)
end

"variable: `vmc[j]` for `j` in `convdc`"
function variable_converter_internal_voltage_magnitude(pm::AbstractPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded::Bool = true)
    PowerModels.var(pm, nw, cnd)[:vmc] = @variable(pm.model,
    [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_vmc",
    lower_bound = PowerModels.ref(pm, nw, :convdc, i, "Vmmin", cnd),
    upper_bound = PowerModels.ref(pm, nw, :convdc, i, "Vmmax", cnd),
    start = PowerModels.ref(pm, nw, :convdc, i, "Vtar", cnd)
    )
end

"variable: `vac[j]` for `j` in `convdc`"
function variable_converter_internal_voltage_angle(pm::AbstractPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded::Bool = true)
    bigM = 2*pi; #
    PowerModels.var(pm, nw, cnd)[:vac] = @variable(pm.model,
    [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_vac",
    lower_bound = -bigM,
    upper_bound =  bigM,
    start = 0
    )
end

"variable: `wrf_ac[j]` and `wif_ac`  for `j` in `convdc`"
function variable_converter_filter_voltage_cross_products(pm::AbstractPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded::Bool = true)
    bigM = 1.2; # only internal converter voltage is strictly regulated
    PowerModels.var(pm, nw, cnd)[:wrf_ac] = @variable(pm.model,
    [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_wrf_ac",
    lower_bound = 0,
    upper_bound = (PowerModels.ref(pm, nw, :convdc, i, "Vmmax", cnd) * bigM)^2,
    start = PowerModels.ref(pm, nw, :convdc, i, "Vtar", cnd)^2
    )
    PowerModels.var(pm, nw, cnd)[:wif_ac] = @variable(pm.model,
    [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_wif_ac",
    lower_bound = -(PowerModels.ref(pm, nw, :convdc, i, "Vmmax", cnd) * bigM)^2,
    upper_bound =  (PowerModels.ref(pm, nw, :convdc, i, "Vmmax", cnd) * bigM)^2,
    start = PowerModels.ref(pm, nw, :convdc, i, "Vtar", cnd)^2
    )
end

"variable: `wf_ac` for `j` in `convdc`"
function variable_converter_filter_voltage_magnitude_sqr(pm::AbstractPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded::Bool = true)
    bigM = 1.2; # only internal converter voltage is strictly regulated
    PowerModels.var(pm, nw, cnd)[:wf_ac] = @variable(pm.model,
    [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_wf_ac",
    lower_bound =  (PowerModels.ref(pm, nw, :convdc, i, "Vmmin", cnd) / bigM)^2,
    upper_bound =  (PowerModels.ref(pm, nw, :convdc, i, "Vmmax", cnd) * bigM)^2,
    start = PowerModels.ref(pm, nw, :convdc, i, "Vtar", cnd)^2
    )
end

"variable: `wrc_ac[j]` and `wic_ac[j]`  for `j` in `convdc`"
function variable_converter_internal_voltage_cross_products(pm::AbstractPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded::Bool = true)
    bigM = 1.2; # only internal converter voltage is strictly regulated
    PowerModels.var(pm, nw, cnd)[:wrc_ac] = @variable(pm.model,
    [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_wrc_ac",
    lower_bound =  0,
    upper_bound =  (PowerModels.ref(pm, nw, :convdc, i, "Vmmax", cnd) * bigM)^2,
    start = PowerModels.ref(pm, nw, :convdc, i, "Vtar", cnd)^2
    )
    PowerModels.var(pm, nw, cnd)[:wic_ac] = @variable(pm.model,
    [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_wic_ac",
    lower_bound =  -(PowerModels.ref(pm, nw, :convdc, i, "Vmmax", cnd) * bigM)^2,
    upper_bound =   (PowerModels.ref(pm, nw, :convdc, i, "Vmmax", cnd) * bigM)^2,
    start = PowerModels.ref(pm, nw, :convdc, i, "Vtar", cnd)^2
    )
end

"variable: `wc_ac[j]` for `j` in `convdc`"
function variable_converter_internal_voltage_magnitude_sqr(pm::AbstractPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded::Bool = true)
    PowerModels.var(pm, nw, cnd)[:wc_ac] = @variable(pm.model,
    [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_wc_ac",
    lower_bound =  PowerModels.ref(pm, nw, :convdc, i, "Vmmin", cnd)^2,
    upper_bound =  PowerModels.ref(pm, nw, :convdc, i, "Vmmax", cnd)^2,
    start = PowerModels.ref(pm, nw, :convdc, i, "Vtar", cnd)^2
    )
end


text:
using Pkg.Artifacts
using polymake_jll

function prepare_deps_tree()
   mutable_artifacts_toml = joinpath(dirname(pathof(polymake_jll)), "..", "MutableArtifacts.toml")
   polymake_tree = "polymake_tree"
   polymake_tree_hash = artifact_hash(polymake_tree, mutable_artifacts_toml)

   # create a directory tree for polymake with links to dependencies
   # looking similiar to the tree in the build environment
   # for compiling wrappers at run-time
   polymake_tree_hash = create_artifact() do art_dir
      mkpath(joinpath(art_dir,"deps"))
      for dep in [polymake_jll.FLINT_jll,
                  polymake_jll.GMP_jll,
                  polymake_jll.MPFR_jll,
                  polymake_jll.PPL_jll,
                  polymake_jll.Perl_jll,
                  polymake_jll.bliss_jll,
                  polymake_jll.boost_jll,
                  polymake_jll.cddlib_jll,
                  polymake_jll.lrslib_jll,
                  polymake_jll.normaliz_jll]
         symlink(dep.artifact_dir, joinpath(art_dir,"deps","$dep"))
      end
      for dir in readdir(polymake_jll.artifact_dir)
         symlink(joinpath(polymake_jll.artifact_dir,dir), joinpath(art_dir,dir))
      end
   end

   bind_artifact!(mutable_artifacts_toml,
      polymake_tree,
      polymake_tree_hash;
      force=true
   )

   # Point polymake to our custom tree
   return artifact_path(polymake_tree_hash)
end

max_stars_repo_name: vchuravy/Yggdrasil

max_issues_repo_name: vchuravy/Yggdrasil

max_forks_repo_name: vchuravy/Yggdrasil


text:
include("args.jl")

function sw_p_to_lambda_den(sw, p)
    sw = tf.reshape(sw, (1, m, n, 1))
    p = tf.reshape(p, (1, m, n, 1))
    sw = tf.image.resize_bilinear(sw, (nz, nx))
    p = tf.image.resize_bilinear(p, (nz, nx))
    sw = cast(sw, Float64)
    p = cast(p, Float64)
    sw = squeeze(sw)
    p = squeeze(p)
    # tran_lambda, tran_den = Gassman(sw)
    # tran_lambda, tran_den = RockLinear(sw) # test linear relationship
    tran_lambda, tran_den = Patchy(sw)
    # tran_lambda_pad =  tf.pad(tran_lambda, [nPml (nPml+nPad); nPml nPml], constant_values=3500.0^2*2200.0/3.0) /1e6
    # tran_den_pad = tf.pad(tran_den, [nPml (nPml+nPad); nPml nPml], constant_values=2200.0)
    return tran_lambda, tran_den
end

args["version"] = "CO2"
lambdasObs = Array{PyObject}(undef, n_survey)
densObs = Array{PyObject}(undef, n_survey)
for iSur = 2:n_survey
  lp = readdlm("./$(args["version"])/FWI_stage$(iSur)/loss.txt")
  Lp = Int64((lp[end,1]))
  lambdasObs[iSur] = constant(readdlm("./$(args["version"])/FWI_stage$(iSur)/Lambda$Lp.txt")*1e6)
  densObs[iSur] = constant(readdlm("./$(args["version"])/FWI_stage$(iSur)/Den$Lp.txt")*1e6)
end

tfCtxInit = tfCtxGen(m,n,h,NT,Δt,Z,X,ρw,ρo,μw,μo,K_init,g,ϕ,qw,qo, sw0, false)
out_sw_init, out_p_init = imseq(tfCtxInit)
lambdas = Array{PyObject}(undef, n_survey)
dens = Array{PyObject}(undef, n_survey)
for i = 1:n_survey
    sw = out_sw_init[survey_indices[i]]
    p = out_p_init[survey_indices[i]]
    lambdas[i], dens[i] = sw_p_to_lambda_den(sw, p)
end
lambdasObs[1] = lambdas[1]
densObs[1] = dens[1]

function objective_function(lambdasObs, lambdas, densObs, dens)
    # tf.nn.l2_loss(lambdasObs - lambdas) + tf.nn.l2_loss(densObs - dens)
    # tf.nn.l2_loss(densObs - dens)
    loss = constant(0.0)
    for i=1:n_survey
        loss += tf.nn.l2_loss(lambdasObs[i][10:nz-10,10:nx-10] - lambdas[i][10:nz-10,10:nx-10])
    end
    return loss
end

J = objective_function(lambdasObs, lambdas, densObs, dens)
gradK = gradients(J, tfCtxInit.K)

if !isdir("./$(args["version"])/flow_fit_results")
    mkdir("./$(args["version"])/flow_fit_results")
end

__cnt = 0
# invK = zeros(m,n)
function print_loss(l, invK, grad)
    global __cnt, __l, __K, __grad
    if mod(__cnt,1)==0
        println("\niter=$__iter, eval=$__cnt, current loss=",l)
        # println("a=$a, b1=$b1, b2=$b2")
    end
    __l = l
    __K = invK
    __grad = grad
    __cnt += 1
end

__iter = 0
function print_iter(rk)
    global __iter, __l, __K, __grad
    if mod(__iter,1)==0
        println("\n************* ITER=$__iter *************\n")
    end
    __iter += 1
    open("./$(args["version"])/flow_fit_results/loss.txt", "a") do io
        writedlm(io, Any[__iter __l])
    end
    writedlm("./$(args["version"])/flow_fit_results/K$__iter.txt", __K)
    writedlm("./$(args["version"])/flow_fit_results/gradK$__iter.txt", __grad)
end
# u = readdlm("...") --> plot(u[:,2])

# config = tf.ConfigProto()
# config.intra_op_parallelism_threads = 24
# config.inter_op_parallelism_threads = 24
sess = Session(); init(sess);
opt = ScipyOptimizerInterface(J, var_list=[tfCtxInit.K], var_to_bounds=Dict(tfCtxInit.K=> (10.0, 130.0)), method="L-BFGS-B",
    options=Dict("maxiter"=> 100, "ftol"=>1e-12, "gtol"=>1e-12))
ScipyOptimizerMinimize(sess, opt, loss_callback=print_loss, step_callback=print_iter, fetches=[J, tfCtxInit.K, gradK])


text:
module ClearMethods

const MODULECOLORS = [
    :light_blue,
    :light_yellow,
    :light_red,
    :light_green,
    :light_magenta,
    :light_cyan,
    :blue,
    :yellow,
    :red,
    :green,
    :magenta,
    :cyan,
]
const EXPAND_BASE_PATHS = Ref(true)
const CONTRACT_USER_DIR = Ref(true)

function getfile(file)
    if EXPAND_BASE_PATHS[]
        file = expandbasepath(file)
    end
    if CONTRACT_USER_DIR[]
        file = replaceuserpath(file)
    end
    return file
end

function expandbasepath(str)
    basefileregex = if Sys.iswindows()
        r"^\.\\\w+\.jl$"
    else
        r"^\./\w+\.jl$"
    end
    if !isnothing(match(basefileregex, str))
        sourcestring = Base.find_source_file(str[3:end]) # cut off ./
    else
        str
    end
end

function replaceuserpath(str)
    str1 = replace(str, homedir() => "~")
    # seems to be necessary for some paths with small letter drive c:// etc
    replace(str1, lowercasefirst(homedir()) => "~")
end

function Base.show(io::IO, m::Method, modulecolor = :yellow)
    tv, decls, file, line = Base.arg_decl_parts(m)
    sig = Base.unwrap_unionall(m.sig)
    if sig === Tuple
        # Builtin
        print(io, m.name, "(...) in ", m.module)
        return
    end
    printstyled(io, m.name; bold = true)
    if !isempty(decls[2:end])
        # type signature
        printstyled(io, "(", color = :light_black)

        for (i, (varname, vartype)) in enumerate(decls[2:end])
            if i > 1
                printstyled(io, ", ", color = :light_black)
            end
            printstyled(io, string(varname), color = :light_black, bold = true)
            if !isempty(vartype)
                printstyled(io, "::")
                printstyled(io, string(vartype), color = :light_black)
            end
        end
    end
    kwargs = Base.kwarg_decl(m)
    if !isempty(kwargs)
        print(io, "; ")
        for k in kwargs
            printstyled(io, string(k), color = :light_black, bold = true)
            print(io, ", ")
        end
    end
    printstyled(io, ")", color = :light_black)
    show_method_params(io, tv)
    println(io)
    printstyled(io, " " * "@ ", color = :light_black)
    # if !isempty(m.module)
    #     printstyled(io, m.module, color = :yellow)
    #     print(io, " ")
    # end
    printstyled(io, m.module, color = modulecolor)
    print(io, " ")
    if line > 0
        file, line = Base.updated_methodloc(m)
        file = getfile(file)
    end

    pathparts = splitpath(file)
    folderparts = pathparts[1:end-1]
    if !isempty(folderparts)
        printstyled(
            io,
            joinpath(folderparts...) * (Sys.iswindows() ? "\\" : "/"),
            color = :light_black,
        )
    end

    # filename, separator, line
    # bright black (90) and underlined (4)
    print(io, "\033[90;4m$(pathparts[end] * ":" * string(line))\033[0m")
end

function show_method_params(io::IO, tv)
    if !isempty(tv)
        printstyled(io, " where "; color = :light_black)
        if length(tv) == 1
            printstyled(io, tv[1]; color = :light_black)
        else
            printstyled(io, "{"; color = :light_black)
            for i = 1:length(tv)
                if i > 1
                    printstyled(io, ", "; color = :light_black)
                end
                x = tv[i]
                printstyled(io, x; color = :light_black)
                io = IOContext(io, :unionall_env => x)
            end
            printstyled(io, "}"; color = :light_black)
        end
    end
end # function show_method_params

getmodule(m::Method) = m.module
getmodule(list::Base.MethodList) = map(getmodule, list)

function Base.show_method_table(
    io::IO,
    ms::Base.MethodList,
    max::Int = -1,
    header::Bool = true,
)
    mt = ms.mt
    name = mt.name
    hasname = isdefined(mt.module, name) && typeof(getfield(mt.module, name)) <: Function
    if header
        Base.show_method_list_header(io, ms, str -> "\"" * str * "\"")
    end
    n = rest = 0
    local last
    LAST_SHOWN_LINE_INFOS = get(io, :LAST_SHOWN_LINE_INFOS, Tuple{String,Int}[])

    modules = getmodule(ms)

    uniquemodules = setdiff(unique(modules), [""])
    modulecolors = Dict(
        u => c for (u, c) in Iterators.zip(uniquemodules, Iterators.cycle(MODULECOLORS))
    )

    resize!(LAST_SHOWN_LINE_INFOS, 0)
    for meth in ms
        if max == -1 || n < max
            n += 1
            println(io)
            print(io, "[$n] ")
            modulecolor = get(modulecolors, meth.module, :default)
            show(io, meth, modulecolor)
            file, line = Base.updated_methodloc(meth)
            push!(LAST_SHOWN_LINE_INFOS, (string(file), line))
        else
            rest += 1
            last = meth
        end
    end
    if rest > 0
        println(io)
        if rest == 1
            show(io, last)
        else
            print(io, "... $rest methods not shown")
            if hasname
                print(io, " (use methods($name) to see them all)")
            end
        end
    end
end

end


text:
#!/usr/bin/env julia
# Copyright (c) 2016, Jost von Hardenberg - ISAC-CNR, Italy

using RainFARM
using ArgParse
using Printf
using DelimitedFiles

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "infile"
            help = "The input file to downscale"
            arg_type = AbstractString
            required = true
        "--varname", "-v"
            help = "Input variable name"
            arg_type = AbstractString
            default = ""
        "--outfile", "-o", "--out"
            help = "Output filename for spectrum"
            arg_type = AbstractString
            default = ""
        "--kmin", "-k"
            help = "Minimum wavenumber"
            arg_type = Int
            default = 1
    end

    s.description = "Estimation of spatial spectral slope for RainFARM downscaling"

    return parse_args(s)
end

args = parse_commandline()
filenc = args["infile"]
varnc = args["varname"]
outfile = args["outfile"]
kmin = args["kmin"]

#println("Estimating slope ",filenc)

(pr, lon_mat, lat_mat) = read_netcdf2d(filenc, varnc)
#println("Size var:", size(pr)," size lon: ",size(lon_mat)," size lat: ", size(lat_mat))
# Compute fft and slope
(fxp, ftp, fs) = fft3d(pr)
sx = fitslope(fxp, kmin=kmin)
#println("Computed spatial spectral slope: ",sx)
println(sx)

if outfile!=""
   nk = length(fxp[:])
   k = collect(1:nk)
   aa = zeros(nk, 2)
   aa[:,1] = k
   aa[:,2] = fxp
   fname = "$outfile.s.dat"
   writedlm(fname, aa)
   fname = "$outfile.t.dat"
   nw = length(ftp[:])
   w = collect(1:nw)
   aa = zeros(nw,2)
   aa[:,1] = w
   aa[:,2] = ftp
   writedlm(fname, aa)
   fname = "$outfile.2d.dat"
   writedlm(fname, fs)
end


text:
# Mel-generalized cepstrum analysis

function mcep!(mc::StridedVector{Cdouble}, windowed::StridedVector{Cdouble},
               α=0.35;
               miniter::Int=2,
               maxiter::Int=30,
               threshold::Float64=0.001,
               etype::Int=0,
               eps::Float64=0.0,
               min_det::Float64=1.0e-6,
               itype::Int=0)
    if itype ∉ 0:4
        throw(ArgumentError("unsupported itype: $itype, must be ∈ 0:4"))
    end
    if etype ∉ 0:2
        throw(ArgumentError("unsupported etype: $etype, must be ∈ 0:2"))
    end
    if etype == 0 && eps != 0.0
        throw(ArgumentError("eps cannot be specified for etype = 0"))
    end
    if etype ∈ 1:2 && eps < 0.0
        throw(ArgumentError("eps: $eps, must be >= 0"))
    end
    if min_det < 0.0
        throw(ArgumentError("min_det must be positive: min_det = $min_det"))
    end

    order = length(mc) - 1
    ret = ccall((:mcep, libSPTK), Cint,
                (Ptr{Cdouble}, Cint, Ptr{Cdouble}, Cint,
                 Cdouble, Cint, Cint, Cdouble, Cint, Cdouble, Cdouble, Cint),
                windowed, length(windowed), mc, order, α,
                miniter, maxiter, threshold, etype, eps, min_det, itype)
    @assert ret ∈ -1:0 || ret ∈ 3:4
    if ret == 3
        error("failed to compute mcep; error occured in theq")
    elseif ret == 4
        error("zero(s) are found in periodogram, use eps option to floor")
    end

    mc
end

function mcep(windowed::StridedVector{Cdouble}, order=25, α=0.35; kargs...)
    mc = Array{Cdouble}(order+1)
    mcep!(mc, windowed, α; kargs...)
end

function gcep!(gc::StridedVector{Cdouble}, windowed::StridedVector{Cdouble},
               γ=0.0;
               miniter::Int=2,
               maxiter::Int=30,
               threshold::Float64=0.001,
               etype::Int=0,
               eps::Float64=0.0,
               min_det::Float64=1.0e-6,
               itype::Int=0,
               norm::Bool=false)
    assert_gamma(γ)
    if itype ∉ 0:4
        throw(ArgumentError("unsupported itype: $itype, must be ∈ 0:4"))
    end
    if etype ∉ 0:2
        throw(ArgumentError("unsupported etype: $etype, must be ∈ 0:2"))
    end
    if etype == 0 && eps != 0.0
        throw(ArgumentError("eps cannot be specified for etype = 0"))
    end
    if etype ∈ 1:2 && eps < 0.0
        throw(ArgumentError("eps: $eps, must be >= 0"))
    end
    if min_det < 0.0
        throw(ArgumentError("min_det must be positive: min_det = $min_det"))
    end

    order = length(gc) - 1
    # Note that the output is gain-normalized
    ret = ccall((:gcep, libSPTK), Cint,
                (Ptr{Cdouble}, Cint, Ptr{Cdouble}, Cint,
                 Cdouble, Cint, Cint, Cdouble, Cint, Cdouble, Cdouble, Cint),
                windowed, length(windowed), gc, order, γ,
                miniter, maxiter, threshold, etype, eps, min_det, itype)
    @assert ret ∈ -1:0 || ret == 3
    if ret == 3
        error("failed to compute gcep; error occured in theq")
    end

    !norm && ignorm!(gc, gc, γ)

    gc
end

function gcep(windowed::StridedVector{Cdouble}, order=25, γ=0.0; kargs...)
    gc = Array{Cdouble}(order + 1)
    gcep!(gc, windowed, γ; kargs...)
end

function mgcep!(mgc::StridedVector{Cdouble}, windowed::StridedVector{Cdouble},
                α=0.35,
                γ=0.0;
                num_recursions::Int=length(windowed)-1,
                miniter::Int=2,
                maxiter::Int=30,
                threshold::Float64=0.001,
                etype::Int=0,
                eps::Float64=0.0,
                min_det::Float64=1.0e-6,
                itype::Int=0,
                otype::Int=0)
    assert_gamma(γ)
    if itype ∉ 0:4
        throw(ArgumentError("unsupported itype: $itype, must be ∈ 0:4"))
    end
    if etype ∉ 0:2
        throw(ArgumentError("unsupported etype: $etype, must be ∈ 0:2"))
    end
    if etype == 0 && eps != 0.0
        throw(ArgumentError("eps cannot be specified for etype = 0"))
    end
    if etype ∈ 1:2 && eps < 0.0
        throw(ArgumentError("eps: $eps, must be >= 0"))
    end
    if min_det < 0.0
        throw(ArgumentError("min_det must be positive: min_det = $min_det"))
    end
    if otype ∉ 0:5
        throw(ArgumentError("unsupported otype: $otype, must be ∈ 0:5"))
    end

    order = length(mgc) - 1
    ret = ccall((:mgcep, libSPTK), Cint,
                (Ptr{Cdouble}, Cint, Ptr{Cdouble}, Cint, Cdouble, Cdouble,
                 Cint, Cint, Cint, Cdouble, Cint, Cdouble, Cdouble, Cint),
                windowed, length(windowed), mgc, order, α, γ, num_recursions,
                miniter, maxiter, threshold, etype, eps, min_det, itype)
    @assert ret ∈ -1:0 || ret == 3
    if ret == 3
        error("failed to compute mgcep; error occured in theq")
    end

    if otype ∈ 0:2 || otype == 4
        ccall((:ignorm, libSPTK), Void, (Ptr{Cdouble}, Ptr{Cdouble},
                                         Cint, Cdouble),
              mgc, mgc, order, γ)
    end

    if otype == 0 || otype == 2 || otype == 4
        ccall((:b2mc, libSPTK), Void,
              (Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cdouble),
              mgc, mgc, order, α)
    end

    if otype == 2 || otype == 4
        ccall((:gnorm, libSPTK), Void, (Ptr{Cdouble}, Ptr{Cdouble},
                                        Cint, Cdouble),
              mgc, mgc, order, γ)
    end

    if otype ∈ 4:5
        for i=2:length(mgc)
            mgc[i] *= γ
        end
    end

    mgc
end

function mgcep(windowed::StridedVector{Cdouble}, order=25, α=0.35, γ=0.0;
               kargs...)
    mgc = Array{Cdouble}(order + 1)
    mgcep!(mgc, windowed, α, γ; kargs...)
end

function uels!(c::StridedVector{Cdouble}, windowed::StridedVector{Cdouble};
               miniter::Int=2,
               maxiter::Int=30,
               threshold::Float64=0.001,
               etype::Int=0,
               eps::Float64=0.0,
               itype::Int=0)
    if itype ∉ 0:4
        throw(ArgumentError("unsupported itype: $itype, must be ∈ 0:4"))
    end
    if etype ∉ 0:2
        throw(ArgumentError("unsupported etype: $etype, must be ∈ 0:2"))
    end
    if etype == 0 && eps != 0.0
        throw(ArgumentError("eps cannot be specified for etype = 0"))
    end
    if etype ∈ 1:2 && eps < 0.0
        throw(ArgumentError("eps: $eps, must be >= 0"))
    end

    order = length(c) - 1
    ret = ccall((:uels, libSPTK), Cint,
                (Ptr{Cdouble}, Cint, Ptr{Cdouble}, Cint, Cint, Cint, Cdouble, Cint,
                 Cdouble, Cint),
                windowed, length(windowed), c, order,
                miniter, maxiter, threshold, etype, eps, itype)
    @assert ret ∈ -1:0 || ret == 3
    if ret == 3
        error("zero(s) are found in periodogram, use eps option to floor")
    end

    c
end

function uels(windowed::StridedVector{Cdouble}, order=25; kargs...)
    c = Array{Cdouble}(order + 1)
    uels!(c, windowed; kargs...)
end

function fftcep!(c::StridedVector{Cdouble}, logsp::StridedVector{Cdouble};
                 num_iter::Int=0,
                 acceleration_factor::Float64=0.0)
    order = length(c) - 1
    ccall((:fftcep, libSPTK), Void,
          (Ptr{Cdouble}, Cint, Ptr{Cdouble}, Cint, Cint, Cdouble),
          logsp, length(logsp), c, order, num_iter, acceleration_factor)
    c
end

function fftcep(logsp::StridedVector{Cdouble}, order=25; kargs...)
    c = Array{Cdouble}(order + 1)
    fftcep!(c, logsp; kargs...)
end

function lpc!(a::StridedVector{Cdouble}, x::StridedVector{Cdouble};
              min_det::Float64=1.0e-6)
    if min_det < 0.0
        throw(ArgumentError("min_det must be positive: min_det = $min_det"))
    end

    order = length(a) - 1
    ret = ccall((:lpc, libSPTK), Cint,
                (Ptr{Cdouble}, Cint, Ptr{Cdouble}, Cint, Cdouble),
                x, length(x), a, order, min_det)
    @assert ret ∈ -2:0
    if ret == -2
        warn("failed to compute `stable` LPC. Please try again with different parameters")
    elseif ret == -1
        error("failed to compute LPC. Please try again with different parameters")
    end

    a
end

function lpc(x::StridedVector{Cdouble}, order=25; kargs...)
    a = Array{Cdouble}(order+1)
    lpc!(a, x; kargs...)
end


text:
"Check that vect is correctly specified"
function chkvect(vect::AbstractChar)
    if !(vect == 'N' || vect == 'V')
        throw(ArgumentError("vect argument must be 'N' (X is not returned), 'V' (X is returned)"))
    end
    vect
end

## Symmetric tridiagonalization
for (fname, elty) in ((:dsbtrd_,:Float64),
                      (:ssbtrd_,:Float32))
    @eval begin
        function sbtrd!(vect::Char, uplo::Char,
                        m::Int, k::Int, A::AbstractMatrix{$elty},
                        d::AbstractVector{$elty}, e::AbstractVector{$elty}, q::AbstractVector{$elty},
                        work::AbstractVector{$elty})
            @assert !has_offset_axes(A)
            chkstride1(A)
            chkuplo(uplo)
            chkvect(vect)
            info = Ref{BlasInt}()
            n    = size(A,2)
            n ≠ m && throw(ArgumentError("Matrix must be square"))
            size(A,1) < k+1 && throw(ArgumentError("Not enough bands"))
            info  = Ref{BlasInt}()
            ccall((@blasfunc($fname), liblapack), Nothing,
                (Ref{UInt8}, Ref{UInt8},
                 Ref{BlasInt}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                 Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ptr{BlasInt}),
                 vect, uplo,
                 n, k, A, max(1,stride(A,2)),
                 d, e, q, max(1,stride(q,2)), work, info)
            LAPACK.chklapackerror(info[])
            d, e, q
        end
    end
end

# (GB) general banded matrices, LU decomposition
for (gbtrf, elty) in
    ((:dgbtrf_, :Float64),
     (:sgbtrf_, :Float32),
     (:zgbtrf_, :ComplexF64),
     (:cgbtrf_, :ComplexF32))
    @eval begin
        # SUBROUTINE DGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO )
        # *     .. Scalar Arguments ..
        #       INTEGER            INFO, KL, KU, LDAB, M, N
        # *     .. Array Arguments ..
        #       INTEGER            IPIV( * )
        #       DOUBLE PRECISION   AB( LDAB, * )
        function gbtrf!(m::Int, l::Int, u::Int, A::AbstractMatrix{$elty},
                        ipiv::AbstractVector{BlasInt})
            @assert !has_offset_axes(A)
            chkstride1(A)
            info = Ref{BlasInt}()
            n    = size(A,2)
            size(A,1) < l+u+1 && throw(ArgumentError("Not enough bands"))
            length(ipiv) < min(m,n) && throw(ArgumentError("Not enough pivots"))
            info = Ref{BlasInt}()
            ccall((@blasfunc($gbtrf), liblapack), Nothing,
                  (Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
                   Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}),
                  m, n, l, u,
                  A, max(1,stride(A,2)), ipiv, info)
            LAPACK.chklapackerror(info[])
            A, ipiv
        end
    end
end

#     gbtrf!(m, kl, ku, AB) -> (AB, ipiv)
# Compute the LU factorization of a banded matrix `AB`. `kl` is the first
# subdiagonal containing a nonzero band, `ku` is the last superdiagonal
# containing one, and `m` is the first dimension of the matrix `AB`. Returns
# the LU factorization in-place and `ipiv`, the vector of pivots used.
# AB should be a (2l+u) by n matrix.

function gbtrf!(m::Int, l::Int, u::Int, AB::AbstractMatrix)
    ipiv = similar(AB, BlasInt, min(m, size(AB,2)))
    gbtrf!(m, l, u, AB, ipiv)
end

# (GB) general banded matrices, LU solver

#     gbtrs!(trans, kl, ku, m, AB, ipiv, B)
# Solve the equation `AB * X = B`. `trans` determines the orientation of `AB`. It may
# be `N` (no transpose), `T` (transpose), or `C` (conjugate transpose). `kl` is the
# first subdiagonal containing a nonzero band, `ku` is the last superdiagonal
# containing one, and `m` is the first dimension of the matrix `AB`. `ipiv` is the vector
# of pivots returned from `gbtrf!`. Returns the vector or matrix `X`, overwriting `B` in-place.

for (gbtrs, elty) in
    ((:dgbtrs_,:Float64),
     (:sgbtrs_,:Float32),
     (:zgbtrs_,:ComplexF64),
     (:cgbtrs_,:ComplexF32))
    @eval begin

        # SUBROUTINE DGBTRS( TRANS, N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB, INFO)
        # *     .. Scalar Arguments ..
        #       CHARACTER          TRANS
        #       INTEGER            INFO, KL, KU, LDAB, LDB, N, NRHS
        # *     .. Array Arguments ..
        #       INTEGER            IPIV( * )
        #       DOUBLE PRECISION   AB( LDAB, * ), B( LDB, * )
        function gbtrs!(trans::Char, m::Int, l::Int, u::Int,
                        A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt},
                        B::AbstractVecOrMat{$elty})
            @assert !has_offset_axes(A)
            chkstride1(A)
            chktrans(trans)
            info = Ref{BlasInt}()
            n    = size(A,2)
            if m != n || m != size(B,1)
                throw(DimensionMismatch("matrix A has dimensions $(size(A)), but right hand side matrix B has dimensions $(size(B))"))
            end
            size(A,1) < 2*l+u+1 && throw(ArgumentError("not enough bands"))
            info = Ref{BlasInt}()

            ccall((@blasfunc($gbtrs), liblapack), Nothing,
                  (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
                   Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt},
                   Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}),
                  trans, m, l, u, size(B,2),
                  A, max(1,stride(A,2)), ipiv,
                  B, max(1,stride(B,2)), info)
            LAPACK.chklapackerror(info[])
            B
        end
    end
end

# Symmetric/Hermitian Positive Definite banded Cholesky factorization
for (fname, elty) in ((:dpbtrf_,:Float64),
                      (:spbtrf_,:Float32),
                      (:zpbtrf_,:ComplexF64),
                      (:cpbtrf_,:ComplexF32))
    @eval begin
                # SUBROUTINE DPBTRF( UPLO, N, KD, AB, LDAB, INFO )
                # CHARACTER          UPLO
                # INTEGER            INFO, KD, LDAB, N
                # DOUBLE PRECISION   AB( LDAB, * )

        function pbtrf!(uplo::Char, m::Int, kd::Int, A::AbstractMatrix{$elty})
            @assert !has_offset_axes(A)
            chkstride1(A)
            chkuplo(uplo)
            info = Ref{BlasInt}()
            n    = size(A,2)
            n ≠ m && throw(ArgumentError("Matrix must be square"))
            size(A,1) < kd+1 && throw(ArgumentError("Not enough bands"))
            ccall((@blasfunc($fname), liblapack), Nothing,
                (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                 Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}),
                 uplo, n, kd,
                 A, max(1,stride(A,2)), info)
            LAPACK.chklapackerror(info[])
            A
        end
    end
end

# Symmetric/Hermitian Positive Definite split-Cholesky factorization
for (fname, elty) in ((:dpbstf_,:Float64),
                      (:spbstf_,:Float32),
                      (:zpbstf_,:ComplexF64),
                      (:cpbstf_,:ComplexF32))
    @eval begin
                # SUBROUTINE DPBSTF( UPLO, N, KD, AB, LDAB, INFO )
                # CHARACTER          UPLO
                # INTEGER            INFO, KD, LDAB, N
                # DOUBLE PRECISION   AB( LDAB, * )

        function pbstf!(uplo::Char, m::Int, kd::Int, A::AbstractMatrix{$elty})
            @assert !has_offset_axes(A)
            chkstride1(A)
            chkuplo(uplo)
            info = Ref{BlasInt}()
            n    = size(A,2)
            n ≠ m && throw(ArgumentError("Matrix must be square"))
            size(A,1) < kd+1 && throw(ArgumentError("Not enough bands"))
            ccall((@blasfunc($fname), liblapack), Nothing,
                (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                 Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}),
                 uplo, n, kd,
                 A, max(1,stride(A,2)), info)
            LAPACK.chklapackerror(info[])
            A
        end
    end
end

# Convert Symmetric Positive Definite generalized eigenvalue problem
# to a symmetric eigenvalue problem assuming B has been processed by
# a split-Cholesky factorization.
for (fname, elty) in ((:dsbgst_,:Float64),
                      (:ssbgst_,:Float32))
    @eval begin
                # SUBROUTINE DSBGST( VECT, UPLO, N, KA, KB, AB, LDAB,
                # BB, LDBB, X, LDX, WORK, INFO )
                # CHARACTER          UPLO, VECT
                # INTEGER            INFO, KA, KB, LDAB, LDBB, LDX, N
                # DOUBLE PRECISION   AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
                # X( LDX, * )

        function sbgst!(vect::Char, uplo::Char, n::Int, ka::Int, kb::Int,
                         AB::AbstractMatrix{$elty}, BB::AbstractMatrix{$elty},
                         X::AbstractVecOrMat{$elty}, work::AbstractVector{$elty})
            @assert !has_offset_axes(AB, BB, X, work)
            chkstride1(AB, BB)
            chkuplo(uplo)
            chkvect(vect)
            info  = Ref{BlasInt}()
            ccall((@blasfunc($fname), liblapack), Nothing,
                (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                 Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty},
                 Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty},
                 Ptr{BlasInt}),
                 vect, uplo, n, ka,
                 kb, AB, max(stride(AB,2),1), BB,
                 max(stride(BB,2),1), X, max(1,stride(X,2)), work,
                 info)
            LAPACK.chklapackerror(info[])
            AB
        end
    end
end


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

struct URLParsingError <: Exception
    msg::String
end
Base.show(io::IO, p::URLParsingError) = println("HTTP.URLParsingError: ", p.msg)

struct Offset
    off::UInt16
    len::UInt16
end
Offset() = Offset(0, 0)
Base.getindex(A::Vector{UInt8}, o::Offset) = A[o.off:(o.off + o.len - 1)]
Base.isempty(o::Offset) = o.off == 0x0000 && o.len == 0x0000
==(a::Offset, b::Offset) = a.off == b.off && a.len == b.len
const EMPTYOFFSET = Offset()

@enum(http_parser_url_fields,
      UF_SCHEME   = 1
    , UF_HOSTNAME = 2
    , UF_PORT     = 3
    , UF_PATH     = 4
    , UF_QUERY    = 5
    , UF_FRAGMENT = 6
    , UF_USERINFO = 7
    , UF_MAX      = 8
)
const UF_SCHEME_MASK = 0x01
const UF_HOSTNAME_MASK = 0x02
const UF_PORT_MASK = 0x04
const UF_PATH_MASK = 0x08
const UF_QUERY_MASK = 0x10
const UF_FRAGMENT_MASK = 0x20
const UF_USERINFO_MASK = 0x40

@inline function Base.getindex(A::Vector{T}, i::http_parser_url_fields) where {T}
    @inbounds v = A[Int(i)]
    return v
end
@inline function Base.setindex!(A::Vector{T}, v::T, i::http_parser_url_fields) where {T}
    @inbounds v = setindex!(A, v, Int(i))
    return v
end

# url parsing
function parseurlchar(s, ch::Char, strict::Bool)
    @anyeq(ch, ' ', '\r', '\n') && return s_dead
    strict && (ch == '\t' || ch == '\f') && return s_dead

    if s == s_req_spaces_before_url
        (ch == '/' || ch == '*') && return s_req_path
        isalpha(ch) && return s_req_schema
    elseif s == s_req_schema
        isalphanum(ch) && return s
        ch == ':' && return s_req_schema_slash
    elseif s == s_req_schema_slash
        ch == '/' && return s_req_schema_slash_slash
        isurlchar(ch) && return s_req_path
    elseif s == s_req_schema_slash_slash
        ch == '/' && return s_req_server_start
        isurlchar(ch) && return s_req_path
    elseif s == s_req_server_with_at
        ch == '@' && return s_dead
        ch == '/' && return s_req_path
        ch == '?' && return s_req_query_string_start
        (isuserinfochar(ch) || ch == '[' || ch == ']') && return s_req_server
    elseif s == s_req_server_start || s == s_req_server
        ch == '/' && return s_req_path
        ch == '?' && return s_req_query_string_start
        ch == '@' && return s_req_server_with_at
        (isuserinfochar(ch) || ch == '[' || ch == ']') && return s_req_server
    elseif s == s_req_path
        (isurlchar(ch) || ch == '@') && return s
        ch == '?' && return s_req_query_string_start
        ch == '#' && return s_req_fragment_start
    elseif s == s_req_query_string_start || s == s_req_query_string
        isurlchar(ch) && return s_req_query_string
        ch == '?' && return s_req_query_string
        ch == '#' && return s_req_fragment_start
    elseif s == s_req_fragment_start
        isurlchar(ch) && return s_req_fragment
        ch == '?' && return s_req_fragment
        ch == '#' && return s
    elseif s == s_req_fragment
        isurlchar(ch) && return s
        (ch == '?' || ch == '#') && return s
    end
    #= We should never fall out of the switch above unless there's an error =#
    return s_dead;
end

function http_parse_host_char(s::http_host_state, ch)
    if s == s_http_userinfo || s == s_http_userinfo_start
        ch == '@' && return s_http_host_start
        isuserinfochar(ch) && return s_http_userinfo
    elseif s == s_http_host_start
        ch == '[' && return s_http_host_v6_start
        ishostchar(ch) && return s_http_host
    elseif s == s_http_host
        ishostchar(ch) && return s_http_host
        ch == ':' && return s_http_host_port_start
    elseif s == s_http_host_v6_end
        ch == ':' && return s_http_host_port_start
    elseif s == s_http_host_v6
        ch == ']' && return s_http_host_v6_end
        (ishex(ch) || ch == ':' || ch == '.') && return s_http_host_v6
        s == s_http_host_v6 && ch == '%' && return s_http_host_v6_zone_start
    elseif s == s_http_host_v6_start
        (ishex(ch) || ch == ':' || ch == '.') && return s_http_host_v6
        s == s_http_host_v6 && ch == '%' && return s_http_host_v6_zone_start
    elseif s == s_http_host_v6_zone
        ch == ']' && return s_http_host_v6_end
        (isalphanum(ch) || @anyeq(ch, '%', '.', '-', '_', '~')) && return s_http_host_v6_zone
    elseif s == s_http_host_v6_zone_start
        (isalphanum(ch) || @anyeq(ch, '%', '.', '-', '_', '~')) && return s_http_host_v6_zone
    elseif s == s_http_host_port || s == s_http_host_port_start
        isnum(ch) && return s_http_host_port
    end
    return s_http_host_dead
end

function http_parse_host(buf, host::Offset, foundat)
    portoff = portlen = uioff = uilen = UInt16(0)
    off = len = UInt16(0)
    s = ifelse(foundat, s_http_userinfo_start, s_http_host_start)

    for i = host.off:(host.off + host.len - 0x0001)
        p = Char(buf[i])
        new_s = http_parse_host_char(s, p)
        new_s == s_http_host_dead && throw(URLParsingError("encountered invalid host character: \n$(String(buf))\n$(lpad("", i-1, "-"))^"))
        if new_s == s_http_host
            if s != s_http_host
                off = i
            end
            len += 0x0001

        elseif new_s == s_http_host_v6
            if s != s_http_host_v6
                off = i
            end
            len += 0x0001

        elseif new_s == s_http_host_v6_zone_start || new_s == s_http_host_v6_zone
            len += 0x0001

        elseif new_s == s_http_host_port
            if s != s_http_host_port
                portoff = i
                portlen = 0x0000
            end
            portlen += 0x0001

        elseif new_s == s_http_userinfo
            if s != s_http_userinfo
                uioff = i
                uilen = 0x0000
            end
            uilen += 0x0001
        end
        s = new_s
    end
    if @anyeq(s, s_http_host_start, s_http_host_v6_start, s_http_host_v6, s_http_host_v6_zone_start,
             s_http_host_v6_zone, s_http_host_port_start, s_http_userinfo, s_http_userinfo_start)
        throw(URLParsingError("ended in unexpected parsing state: $s"))
    end
    # (host, port, userinfo)
    return Offset(off, len), Offset(portoff, portlen), Offset(uioff, uilen)
end

function http_parser_parse_url(buf, startind=1, buflen=length(buf), isconnect::Bool=false)
    s = ifelse(isconnect, s_req_server_start, s_req_spaces_before_url)
    old_uf = UF_MAX
    off = len = 0
    foundat = false
    offsets = Offset[Offset(), Offset(), Offset(), Offset(), Offset(), Offset(), Offset()]
    mask = 0x00
    for i = startind:(startind + buflen - 1)
        @inbounds p = Char(buf[i])
        olds = s
        s = parseurlchar(s, p, false)
        if s == s_dead
            throw(URLParsingError("encountered invalid url character for parsing state = $(ParsingStateCode(olds)): \n$(String(buf))\n$(lpad("", i-1, "-"))^"))
        elseif @anyeq(s, s_req_schema_slash, s_req_schema_slash_slash, s_req_server_start, s_req_query_string_start, s_req_fragment_start)
            continue
        elseif s == s_req_schema
            uf = UF_SCHEME
            mask |= UF_SCHEME_MASK
        elseif s == s_req_server_with_at
            foundat = true
            uf = UF_HOSTNAME
            mask |= UF_HOSTNAME_MASK
        elseif s == s_req_server
            uf = UF_HOSTNAME
            mask |= UF_HOSTNAME_MASK
        elseif s == s_req_path
            uf = UF_PATH
            mask |= UF_PATH_MASK
        elseif s == s_req_query_string
            uf = UF_QUERY
            mask |= UF_QUERY_MASK
        elseif s == s_req_fragment
            uf = UF_FRAGMENT
            mask |= UF_FRAGMENT_MASK
        else
            throw(URLParsingError("ended in unexpected parsing state: $s"))
        end
        if uf == old_uf
            len += 1
            continue
        end
        if old_uf != UF_MAX
            offsets[old_uf] = Offset(off, len)
        end
        off = i
        len = 1
        old_uf = uf
    end
    if old_uf != UF_MAX
        offsets[old_uf] = Offset(off, len)
    end
    check = ~(UF_HOSTNAME_MASK | UF_PATH_MASK)
    if (mask & UF_SCHEME_MASK > 0) && (mask | check == check)
        throw(URLParsingError("URI must include host or path with scheme"))
    end
    if mask & UF_HOSTNAME_MASK > 0
        host, port, userinfo = http_parse_host(buf, offsets[UF_HOSTNAME], foundat)
        if !isempty(host)
            offsets[UF_HOSTNAME] = host
            mask |= UF_HOSTNAME_MASK
        end
        if !isempty(port)
            offsets[UF_PORT] = port
            mask |= UF_PORT_MASK
        end
        if !isempty(userinfo)
            offsets[UF_USERINFO] = userinfo
            mask |= UF_USERINFO_MASK
        end
    end
    # CONNECT requests can only contain "hostname:port"
    if isconnect
        chk = UF_HOSTNAME_MASK | UF_PORT_MASK
        ((mask | chk) > chk) && throw(URLParsingError("connect requests must contain and can only contain both hostname and port"))
    end
    return URI(buf, (offsets[UF_SCHEME],
                     offsets[UF_HOSTNAME],
                     offsets[UF_PORT],
                     offsets[UF_PATH],
                     offsets[UF_QUERY],
                     offsets[UF_FRAGMENT],
                     offsets[UF_USERINFO]))
end


text:
"""
    TransferOperatorTriangulationApprox

A transfer operator estimated over a triangulation of a set of points. Transition
probabilities are computed from approximate simplex intersections.

The return type for the `transferoperator_triangulation_approx` estimator.

## Fields
-  **`transfermatrix`**. The estimated transfer matrix.

"""
struct TransferOperatorTriangulationApprox <: AbstractTransferOperator
    transfermatrix::AbstractArray{Float64, 2}
end

export TransferOperatorTriangulationApprox