text:
export GAP, GAPA, GAPP,  DR, AP, Dykstra, FISTA

include("gap.jl")
include("gapa.jl")
include("gapproj.jl")
include("dykstra.jl")
include("fista.jl")

#Some shortnames for common algorithms
DR(α=0.5; kwargs...) = GAP(α, 2.0, 2.0; kwargs...)
AP(α=1; kwargs...)   = GAP(α, 1.0, 1.0; kwargs...)


text:
#!/usr/bin/env julia

module DBFs

const DBF = Vector{Float32}

export DBF

using Base.Cartesian

"""
use segmentation to get binary image to save memory usage
"""
function compute_DBF(seg::Array{T,3}, obj_id::T) where T
    error("unimplemented")
end

"""

    compute_DBF( pointCloud )

  Returns an array of DBF values for the point cloud. Currently creates
  a binary image, and runs bwd2 on it, though ideally we'd get rid of the
  need for an explicit bin_im
"""
function compute_DBF( pointCloud::Array{T, 2} ) where T
    bin_im = create_binary_image( pointCloud );
    compute_DBF(pointCloud, bin_im)
end

"""
    compute_DBF( bin_im )
"""
function compute_DBF( pointCloud::Array{T,2}, bin_im::Union{BitArray, Array{Bool, 3}} ) where T
    dbf_im = distance_transform( bin_im );
    return extract_dbf_values( dbf_im, pointCloud );
end

"""
compute Distance from Boundary Field (DBF) based on point cloud and the boundary points

WARN: this function do not work correctly!
"""
function compute_DBF( points::Array{T,2}, boundary_point_indexes::Vector ) where T
    error("this function do not work correctly, have a bug!")
    num = size(points, 1)
    dbf = DBF(num)
    fill!(dbf, Inf32)
    for i in 1:num
        point = points[i,:]
        for bpi in boundary_point_indexes
            boundary = points[bpi, :]
            # filter out some far away boundary points
            if  abs(point[1]-boundary[1]) < MAX_BOUNDARY_DISTANCE &&
                abs(point[2]-boundary[2]) < MAX_BOUNDARY_DISTANCE &&
                abs(point[3]-boundary[3]) < MAX_BOUNDARY_DISTANCE
                # compute euclidean distance
                ed = norm(point .- boundary)
                if ed < dbf[i]
                    dbf[i] = ed
                end
            end
        end
    end
    return dbf
end

#NOTE voxelSize not currently functional
"""

    distance_transform( d::AbstractArray{T,N}, voxelSize::Vector{Float32}=ones(Float32, N) )

  Returns a euclidean distance transformation of the mask provided by d. The return
  value will be a volume of the same size as d where the value at each index corresponds
  to the distance between that location and the nearest location for which d > 0.
"""
@generated function distance_transform( d::AbstractArray{T,N}, voxelSize::Vector{Float32}=ones(Float32, N) ) where {T,N}
  quote

  @assert length(voxelSize) == $N;

  res = zeros(Float32,size(d));

  fill_f0!(res, d);

  dims = 1:$N;

  for d in dims
    vol_voronoi_edt!( res, voxelSize[d] );
    res = permutedims( res, circshift(dims,1) );
  end

  sqrt!(res)
  res

  end#quote
end

"""
Fills an n-dimensional volume with initial states for edt transformation,
inf for non-feature voxels, and 0 for feature voxels
"""
@generated function fill_f0!( arr::Array{Float32,N}, fv::AbstractArray{T,N} ) where {T,N}
  quote

  #apparently this generates lots of allocations,
  # I'm still not entirely sure why... but cool!
  # (learned the answer - choice between 0 and Inf is type unstable)
  @nloops $N i arr begin
    (@nref $N arr i) = (@nref $N fv i) > 0 ? 0. : Inf;
  end

  #arr[fv] = 0;
  #arr[!fv] = Inf;

  end#quote
end

"""
Performs the edt transformation along the first dimension of the N-dimensional
volume
"""
@generated function vol_voronoi_edt!( arr::Array{Float32,N}, dim::Float32 ) where N
  quote

  s1 = size(arr,1);
  g = zeros(Float32,(s1,));
  h = zeros(Int,    (s1,));
  @nloops $N i j->(j==1 ? 0 : 1:size(arr,j)) begin

      fill!(h, zero(eltype(h)));
      fill!(g, zero(eltype(g)));
    later_indices = (@ntuple $N j->i_j);
    row_voronoi_edt!( arr, later_indices[2:end], g,h, dim ) #F_{d-1}

  end

  end#quote
end

"""
Performs the edt over a specific row in the volume, following the first dimension
"""
@generated function row_voronoi_edt!( F::Array{Float32,N}, indices::Tuple,
  g::Vector{Float32}, h::Vector{Int}, dim::Float32 ) where N
  quote

  #count of potential feature vectors
  numPotentialFeatureVectors::Int = 0;

  #selecting out the value in the row
  @inbounds f = @nref $N F j->(j==1 ? (:) : indices[j-1]);

  #construct set of feature voxels whose voronoi
  # cells intersect the row
  for i in eachindex( f )

    #scanning for possible feature vector locations
    if !isinf( f[i] )
      if numPotentialFeatureVectors < 2
        #l += 1; g[l] = f[i]; h[l] = i;
        numPotentialFeatureVectors += 1;
        g[numPotentialFeatureVectors] = f[i];
        h[numPotentialFeatureVectors] = i;
      else
        while (numPotentialFeatureVectors >= 2 &&
               remove_euclidean_distance_transform(g[numPotentialFeatureVectors-1],
                                                   g[numPotentialFeatureVectors], f[i],
                                                   h[numPotentialFeatureVectors-1],
                                                   h[numPotentialFeatureVectors], i))
          numPotentialFeatureVectors -= 1;
        end
        numPotentialFeatureVectors += 1;
        g[numPotentialFeatureVectors] = f[i];
        h[numPotentialFeatureVectors] = i;
      end #if numPotentialFeatureVectors
    end #if !isinf

  end #for i

  # if no possible feature voxels, stop now
  if numPotentialFeatureVectors == 0; return nothing end

  #assign new closest feature vectors
  # and update new distances
  num_fvs = numPotentialFeatureVectors; numPotentialFeatureVectors = 1;
  for i in eachindex( f )

    #if we haven't reached the end, and the next feature vector is closer
    # to this location than the current one
    while (numPotentialFeatureVectors < num_fvs &&
           (g[numPotentialFeatureVectors] + (h[numPotentialFeatureVectors] - i)^2 >
                g[numPotentialFeatureVectors+1] + (h[numPotentialFeatureVectors+1] - i)^2))
      numPotentialFeatureVectors += 1;
    end

    (@nref $N F j->(j==1 ? i : indices[j-1])) =  g[numPotentialFeatureVectors] +
                                            (h[numPotentialFeatureVectors] - i)^2;
  end #for i

  end #quote
end

"""
Getting too tired to document these next few, but will be worth it if it works
"""
function fv_isfurther(g1::Float32 ,h1::Int, g2::Float32, h2::Int, i::Int)
  return g1 + (h1 - i)^2 > g2 + (h2 - i)^2;
end

"""
    remove_euclidean_distance_transform
"""
function remove_euclidean_distance_transform( g1::Float32, g2::Float32, g3::Float32, h1::Int, h2::Int, h3::Int )
  a = h2 - h1;
  b = h3 - h2;
  c = h3 - h1;

  return (c*g2 - b*g1 - a*g3 - a*b*c > 0)
end

function sqrt!( d::AbstractArray )
  for i in eachindex(d)
    d[i] = sqrt(d[i]);
  end
end

"""

    create_binary_image( pointCloud )

  Creates a boolean volume where the non-segment indices
  map to true, while the segment indices map to false.
"""
function create_binary_image( pointCloud::Array{T,2} ) where T;

  max_dims = maximum( pointCloud, dims=1 );
    bin_im = falses(max_dims...)

  for p in 1:size( pointCloud, 1 )
    bin_im[ pointCloud[p,:]... ] = false;
  end

  bin_im;
end

"""
    create_binary_image( seg, obj_id )

Creates a boolean volume where the non-segment indices
map to true, while the segment indices map to false
"""
function create_binary_image( seg::Array{T,3}; obj_id::T = one(T) ) where T
    bin_im = trues(size(seg))
    for i in eachindex(seg)
        if seg[i] == obj_id
            bin_im[i] = false
        end
    end
    bin_im
end

"""

    extract_dbf_values( dbf_image, pointCloud )

  Takes an array where rows indicate subscripts, and extracts the values
  within a volume at those subscripts (in row order)
"""
@generated function extract_dbf_values( dbf_image::Array{Float32,N}, pointCloud ) where N
  quote

  num_points = size( pointCloud, 1 );
  dbf_values = zeros(Float32, num_points);

  for p in 1:num_points
    dbf_values[p] = (@nref $N dbf_image i->pointCloud[p,i]);
  end

  dbf_values
  end#quote
end

end#module

ext: jl


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule TracyProfiler_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("TracyProfiler")
JLLWrappers.@generate_main_file("TracyProfiler", UUID("0c351ed6-8a68-550e-8b79-de6f926da83c"))
end  # module TracyProfiler_jll


text:
# When at a leaf node, this function returns the values in an appropriate
# form. This differs for cases when the leaf node holds an array of elements,
# a single strings, or a single number (arrays of numbers not yet implemented).

function get_values(obj)
    if length(obj) > 1
        res = string(obj)[4:end]            # first 3 chars in string are "Any"
    elseif length(obj) == 1
        if typeof(obj[1]) <: AbstractString
            res = string("\"", obj[1], "\"")
        elseif typeof(obj[1]) <: Number
            res = obj[1]
        end
    end
    return res
end

# This function is the workhorse of xml2json(). It uses recursion to fill a
# string with the keys and values of a nested multi-dictionary object
# returned by the xml2dict function.

function tags_and_elements(mdict::DataStructures.MultiDict{Any, Any}, res::ASCIIString, nspaces::Int = 4)
    wspace = get_indentation(nspaces)
    nkeys = length(keys(mdict))
    res *= "\n" * wspace * "{"

    # loop over keys of multi-dictionary, and recursively fill `res`
    # string with keys and values from the multi-dictionary.
    for (i, k) in enumerate(keys(mdict))
        num_elem = length(mdict[k])

        if is_leafnode(mdict[k])
            res *= string("\n", wspace, "\"", k, "\":", get_values(mdict[k]))
            res *= (i < nkeys) ? "," : ""           # comma if non-last entry
        elseif num_elem > 1
            res *= "\n" * wspace * "\"" * k * "\":" * "["
            for j = 1:num_elem
                res = tags_and_elements(mdict[k][j], res, nspaces + 4)
                res *= (j < num_elem) ? "," : ""
            end
            res *= (i < nkeys) ? "\n" * wspace * "]," : "\n" * wspace * "]"
        elseif num_elem == 1
            res *= string("\n", wspace, "\"", k, "\":")
            res = tags_and_elements(mdict[k][1], res, nspaces + 4)
            res *= (i < nkeys) ? "," : ""
        end
    end
    res *= string("\n", wspace, "}")
    return res
end

function xml2json(xroot::LightXML.XMLElement, nspaces::Int = 4, replace_newline::Union{Void, ASCIIString} = "\\n")
    xmlstring = string(xroot)

    # The first element can have attributes, so its name will end either
    # with whitespace or with an close angle bracket (when no attribute).
    angle_bracket_idx = search(xmlstring, ">")[1]
    space_idx = search(xmlstring, " ")[1]
    end_idx = (angle_bracket_idx < space_idx) ? angle_bracket_idx - 1 : space_idx - 1
    heading = string("{", "\n", "\"", xmlstring[2:end_idx], "\"", ":", "\n")

    xdict = xml2dict(xroot, replace_newline)
    body = tags_and_elements(xdict, "", nspaces)
    out_string = heading * body * "\n}"

    return out_string
end


text:
struct VulkanWrapper
    handles::Vector{Expr}
    structs::Vector{Expr}
    funcs::Vector{Expr}
    enums::Vector{Expr}
    docs::Vector{Expr}
end

Base.show(io::IO, vw::VulkanWrapper) = print(
    io,
    "VulkanWrapper with $(length(vw.handles)) handles, $(length(vw.structs)) structs and $(length(vw.funcs)) functions.",
)

include("wrap/classification.jl")
include("wrap/identifiers.jl")
include("wrap/call.jl")
include("wrap/return.jl")
include("wrap/decl.jl")
include("wrap/initialization.jl")
include("wrap/parent.jl")

include("wrap/bitmasks.jl")
include("wrap/overloads.jl")
include("wrap/structs.jl")
include("wrap/functions.jl")
include("wrap/handles.jl")

function VulkanWrapper()
    handles = to_expr.(wrap.(spec_handles))
    structs = to_expr.(wrap.(spec_structs))
    returnedonly_structs = filter(x -> x.is_returnedonly, spec_structs)
    api_structs = filter(!in(returnedonly_structs), spec_structs)
    funcs =
        to_expr.(
            vcat(
                wrap.(spec_funcs),
                add_constructor.(api_structs),
                extend_from_vk.(returnedonly_structs),
                add_constructors.(spec_handles_with_wrappable_constructors)...,
                extend_handle_constructor.(
                    filter(
                        x -> !isnothing(x.create_info_param) && !x.batch,
                        collect(Iterators.flatten(create_funcs.(spec_handles_with_wrappable_constructors))),
                    ),
                ),
                wrap.(spec_funcs; with_func_ptr = true),
                add_constructors.(spec_handles_with_wrappable_constructors; with_func_ptr = true)...,
                extend_handle_constructor.(
                    filter(
                        x -> !isnothing(x.create_info_param) && !x.batch,
                        collect(Iterators.flatten(create_funcs.(spec_handles_with_wrappable_constructors))),
                    );
                    with_func_ptr = true,
                ),
            ),
        )

    enums = to_expr.(wrap.(spec_bitmasks))
    docs = vcat(
        document.(spec_funcs, wrap.(spec_funcs)),
        document.(api_structs, add_constructor.(api_structs)),
        document.(spec_handles_with_wrappable_constructors)...,
    )
    VulkanWrapper(handles, structs, funcs, enums, docs)
end


text:
using NewtonsMethod
using Test

using LinearAlgebra, Statistics, Compat, ForwardDiff

@testset "NewtonsMethod.jl" begin

#This function has no root, it should return nothing
f1(x) = x^2+1
f1_prime(x) = 2*x

@test newtonroot(f1, f1_prime,x₀ = 0).root == nothing
@test newtonroot(f1, x₀ = 0).root == nothing

#Now, I will test some functions which I know the roots
f2(x) = (x - 1)^3
f2_prime(x) = 3*((x-1)^2)

f3(x) = 4*x+5
f3_prime(x) = 4

f4(x) = 5^x - 25
f4_prime(x) =  5^x * log(5)

@test newtonroot(f2, f2_prime,x₀ = 0).root ≈ 0.9999998643434097
@test newtonroot(f2, x₀ = 0).root ≈ 0.9999998643434097

@test newtonroot(f3, f3_prime,x₀ = 0).root ≈ -1.25
@test newtonroot(f3, x₀ = 0).root ≈ -1.25

@test newtonroot(f4, f4_prime,x₀ = 0).root ≈ 2
@test newtonroot(f4, x₀ = 0).root ≈ 2

 #Testing tolerance
 @test !(newtonroot(f2, f2_prime, x₀ = 0, tolerance = 1).root ≈ 0.9999998643434097)
 @test !(newtonroot(f2, x₀ = 0, tolerance = 1).root ≈ 0.9999998643434097)
 @test newtonroot(f2, f2_prime, x₀ = 0, tolerance = 1E-10).root ≈ 1 #with smaller tolerance, I can reach the "exact" number

 #Testing BigFloat
 @test newtonroot(f2, f2_prime, x₀ = BigFloat(0.0), tolerance = 1E-40).root ≈ 1
 @test newtonroot(f2, x₀ = BigFloat(0.0), tolerance = 1E-40).root ≈ 1

 #Testing if the maxiter works as intended
   #Function 2
   @test newtonroot(f2, f2_prime, x₀ = 0, maxiter = 10).root == nothing
  @test newtonroot(f2, x₀ = 0, maxiter = 10).root == nothing
  # Function 4
   @test newtonroot(f4, f4_prime, x₀ = 0, maxiter = 10).root == nothing
   @test newtonroot(f4,x₀ = 0, maxiter = 10).root == nothing

end


text:
#=
Popis
  Funkcia navrati datovu strukturu vrchola zo zadaneho grafu na zaklade jeho indexu.

Parametre
  index : parameter udavajuci index vrchola :: Int64
  g : datova struktura grafu :: GenericGraph

Navratova hodnota
  vert : datova struktura vrchola :: ExVertex

Poznamky
  V pripade ak sa dany vrchol v grafe nenachadza, navratova hodnoa funkcie je false(::Bool).
=#

function get_vertex(index::Int64, g::GenericGraph)

  for vert in vertices(g)
    if vert.index == index
      return vert
    end
  end

  return false

end

max_stars_repo_name: mpuk/GraphAnalysis

max_issues_repo_name: mpuk/GraphAnalysis

max_forks_repo_name: mpuk/GraphAnalysis


text:
#!/usr/bin/env Julia

using Distributions

size, subsize, depth =
    parse(Int, ARGS[1]), parse(Int, ARGS[2]), parse(Int, ARGS[3])

oname = ARGS[4]

function generator(μ, subsize, size, depth)
    N = MvNormal(μ, eye(2) * ((depth * 0.5) ^ 4))
    if depth == 1
        for _ in 1:size
            X = rand(N)
            produce([X' N.μ'])
        end
    else
        for seed in 1:subsize
            generator(rand(N), subsize, size, depth-1)
        end
    end
end

p = Task(() -> generator(zeros(2), subsize, size, depth))

Data = vcat([X for X in p]...)

using DataFrames

df = DataFrame(x=Data[:, 1], y=Data[:, 2])

title = "$(oname)_$(size)_$(subsize)_$(depth)"
writetable("$(title).csv", df)

using Plots

df[:group] = (df[:x] .* df[:y])

plt = scatter(df[:x], df[:y], color=df[:group])

png(title)

max_stars_repo_name: probinso/IR-cluster-rank-demo

max_issues_repo_name: probinso/IR-cluster-rank-demo

max_forks_repo_name: probinso/IR-cluster-rank-demo


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _flux2mag(flux::T, zero_point::T, ABwave::T) where {T<:AbstractFloat}
    if isnan(ABwave)
        return -2.5*log10(flux) - zero_point
    else
        return -2.5*log10(flux) - 5 * log10(float(ABwave)) - 2.406
    end
end

"""
    flux2mag(flux[, zero_point, ABwave=number]) -> magnitude

### Purpose ###

Convert from flux expressed in erg/(s cm² Å) to magnitudes.

### Explanation ###

This is the reverse of `mag2flux`.

### Arguments ###

* `flux`: the flux to be converted in magnitude, expressed in
  erg/(s cm² Å).
* `zero_point`: the zero point level of the magnitude.  If not
 supplied then defaults to 21.1 (Code et al 1976).  Ignored if the `ABwave`
 keyword is supplied
* `ABwave` (optional numeric keyword): wavelength in Angstroms.
 If supplied, then returns Oke AB magnitudes (Oke & Gunn 1983, ApJ, 266, 713;
 http://adsabs.harvard.edu/abs/1983ApJ...266..713O).

### Output ###

The magnitude.

If the `ABwave` keyword is set then magnitude is given by the expression

\$\$\\text{ABmag} = -2.5\\log_{10}(f) - 5\\log_{10}(\\text{ABwave}) - 2.406\$\$

Otherwise, magnitude is given by the expression

\$\$\\text{mag} = -2.5\\log_{10}(\\text{flux}) - \\text{zero point}\$\$

### Example ###

```jldoctest
julia> using AstroLib

julia> flux2mag(5.2e-15)
14.609991640913002

julia> flux2mag(5.2e-15, 15)
20.709991640913003

julia> flux2mag(5.2e-15, ABwave=15)
27.423535345634598
```

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
flux2mag(flux::Real, zero_point::Real=21.1; ABwave::Real=NaN) =
    _flux2mag(promote(float(flux), float(zero_point), float(ABwave))...)


text:
module NewtonSolver

using NLsolve
include("SpectralSolver.jl")
using .SpectralSolver

mutable struct ModeParameters
    l::Int64
    m::Int64
    n::Int64
    s::Int64
    a::Float64
    ω::Complex{Float64}
    Alm::Complex{Float64}
    Nmax::Int64
    lmax::Int64
end

ModeParameters(NT::NamedTuple) = ModeParameters(NT[:l],NT[:m],NT[:n],NT[:s],NT[:a],NT[:ω],NT[:Alm],NT[:Nmax],NT[:lmax])

#Parameter Transformation
function ParameterTransformations(l,m,n,s,a,ω,Alm,Nmax,lmax)
    M = 1
    r₊ = M + sqrt(M^2 - a^2)
    r₋ = M - sqrt(M^2 - a^2)

    σ₊ = (2*ω*M*r₊ - m*a)/(r₊-r₋)
    σ₋ = (2*ω*M*r₋ - m*a)/(r₊-r₋)

    ζ₊ = im*ω
    ξ₋ = (-s-(s+2*im*σ₊))/2
    η₋ = -s+im*σ₋

    ζ = ζ₊
    ξ = ξ₋
    η = η₋

    p = (r₊ - r₋)*(ζ/2)
    α = 1 + s + ξ + η - 2*ζ + s*im*ω/ζ
    γ = 1 + s + 2*η
    δ = 1 + s + 2*ξ
    σ = Alm + (a*ω)^2 - 8*ω^2 + p*(2*α + γ - δ) + (1 + s - (γ + δ)/2 )*(s + (γ + δ)/2 )

    D₀ = δ
    D₁ = 4*p - 2*α + γ - δ -2
    D₂ = 2*α - γ + 2
    D₃ = α*(4*p - δ) - σ
    D₄ = α*(α - γ + 1)

    return ((ζ,ξ,η),(p,α,γ,δ,σ),(D₀,D₁,D₂,D₃,D₄))
end

ParameterTransformations(P::ModeParameters) = ParameterTransformations(P.l,P.m,P.n,P.s,P.a,P.ω,P.Alm,P.Nmax,P.lmax)

function ContinuedFraction(α,β,γ,n::Int,N::Int; rN = zero(α(1)))
    num = β(0)
    for i = 1:n
        num = β(i) - α(i-1)*γ(i)/num
    end
    lhs = num

    num = β(N)+α(N)*rN
    for i = (N):-1:(n+1)
       num = β(i-1) - α(i-1)*γ(i)/num
    end
    rhs = num - β(n)

    lhs + rhs
end

function rNCoeffs(D0,D1,D2,D3,D4)
    u1 = sqrt(-D0 - D1 - D2)
    u2 = (1/2)*(-4 - 2*D0 - D1)
    u3 = (8*u1 + 16*D0*u1 + 8*(D0^2)*u1 + 12*D1*u1 +
     8*D0*D1*u1 + (D1^2)*u1 + 8*D2*u1 +
     4*D0*D2*u1 - 4*D3*u1 - 4*D4*u1)/
   (8*(D0 + D1 + D2))
    u4 = (1/2)*(4 + 4*D0 + 2*D0^2 + D1 + D0*D1 - D3)
    u1,u2,u3,u4
end

function CF(P)
    _,_,(D₀,D₁,D₂,D₃,D₄) = ParameterTransformations(P)
    n = P.n; N = P.Nmax;
    αₙ(n) = n^2 + (D₀+1)*n + D₀
    βₙ(n) = -2*n^2 + (D₁+2)*n +D₃
    γₙ(n) = n^2 + (D₂-3)*n  + D₄ - D₂ + 2
    u1,u2,u3,u4 = rNCoeffs(D₀,D₁,D₂,D₃,D₄)
    rN = 1 + u1*N^(-0.5) + u2*N^(-1) + u3*N^(1.5) + u4*N^(-2)
    ContinuedFraction(αₙ,βₙ,γₙ,n,N; rN = rN)
end

function FiniteDifferenceAtParameters(P,ϵs)
    y = CF(P)
    ωᵢ = P.ω
    Almᵢ = P.Alm

    P.ω = ωᵢ + ϵs+im*0
    P.Alm,_,_ = ComputeAₗₘ(P.s,P.m,P.a*(P.ω), P.Alm,P.lmax)
    yʼr = (CF(P) - y)/ϵs

    P.ω = ωᵢ + im*ϵs
    P.Alm,_,_ = ComputeAₗₘ(P.s,P.m,P.a*(P.ω), P.Alm,P.lmax)
    yʼi = (CF(P) - y)/ϵs

    P.ω = ωᵢ
    P.Alm = Almᵢ
    return yʼr,yʼi,y
end

function Newton!(P,ϵs)
    ω₀ = P.ω
    Alm₀ = P.Alm
    function f!(F,x)
        ωᵢ = P.ω
        Almᵢ = P.Alm
        P.ω = x[1]+im*x[2]
        P.Alm,_,_ = ComputeAₗₘ(P.s,P.m,P.a*(P.ω), P.Alm,P.lmax)
        val = CF(P)
        P.ω = ωᵢ
        P.Alm = Almᵢ
        F[1] = real(val)
        F[2] = imag(val)
    end

    function j!(J,x)
        ω = x[1]+im*x[2]
        P.ω = ω
        ydr,ydi,y = FiniteDifferenceAtParameters(P,ϵs)
        J[1,1] = real(ydr)
        J[1,2] = real(ydi)
        J[2,1] = imag(ydr)
        J[2,2] = imag(ydi)
    end

    function fj!(F,J,x)
        ωᵢ = P.ω
        Almᵢ = P.Alm
        P.ω = x[1]+im*x[2]
        P.Alm,_,_ = ComputeAₗₘ(P.s,P.m,P.a*(P.ω), P.Alm,P.lmax)
        val = CF(P)
        F[1] = real(val)
        F[2] = imag(val)

        P.ω = x[1]+im*x[2] + ϵs+im*0
        P.Alm,_,_ = ComputeAₗₘ(P.s,P.m,P.a*(P.ω), P.Alm,P.lmax)
        yʼr = (CF(P) - val)/ϵs

        P.ω = x[1]+im*x[2] + im*ϵs
        P.Alm,_,_ = ComputeAₗₘ(P.s,P.m,P.a*(P.ω), P.Alm,P.lmax)
        yʼi = (CF(P) - val)/ϵs

        P.ω = ωᵢ
        P.Alm = Almᵢ

        J[1,1] = real(yʼr )
        J[1,2] = real(yʼi)
        J[2,1] = imag(yʼr )
        J[2,2] = imag(yʼi)
    end

    val = CF(P)
    df = OnceDifferentiable(f!, j!, fj!,[real(ω₀),imag(ω₀)],[real(val),imag(val)])
    sol = nlsolve(df, [real(ω₀),imag(ω₀)],ftol = 1e-14, method = :newton)

    ωnew = sol.zero
    P.ω = ωnew[1] + im*ωnew[2]
    Almnew,_,_ = ComputeAₗₘ(P.s,P.m,P.a*(P.ω), P.Alm,P.lmax)
    P.Alm = Almnew
    P
end

export ParameterTransformations, ContinuedFraction, CF, FiniteDifferenceAtParameters, Newton!
export ModeParameters

end

max_stars_repo_name: Potatoasad/AttemptAtQNM

max_issues_repo_name: Potatoasad/AttemptAtQNM

max_forks_repo_name: Potatoasad/AttemptAtQNM


text:
# ScaLapack.jl
module ScaLapack
    # preambles
    include("preambles.jl")
    # scalapack wrapper
    include("scalapack_wrappers.jl")

    # blacs wrapper
    module BLACS
        import ..libscalapack, ..ScaInt, ..f_pchar
        include("blacs.jl")
    end

    # pblas wrapper: planning add thus module ...
    # module PBLAS
    #     import ..libscalapack, ..ScaInt, ..f_pstring
    #     using ..BLACS
    #     include("pblas.jl")
    # end

    # Lite ver. of the ScaLapack
    module ScaLapackLite
        import ..libscalapack, ..ScaInt, ..f_pchar
        import ..zero, ..one, ..desc_idx
        using ..MPI, ..ScaLapack, ..BLACS
        # using ..PBLAS
        include("scalapack_lite.jl")
    end

end


text:
__precompile__()

module ProjectEuler

files = readdir(@__DIR__)
sort!(files)

for i in 2:length(files)
    include(files[i])
end

end # module

max_stars_repo_name: halleysfifthinc/ProjectEuler

max_issues_repo_name: halleysfifthinc/ProjectEuler

max_forks_repo_name: halleysfifthinc/ProjectEuler


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# ╔═╡ 6764c13e-681b-11eb-3527-77caf9d74f5b
using PlutoUI

# ╔═╡ 63a5111a-6817-11eb-2151-fb825d7879ba
using Test

# ╔═╡ d24d5ea0-6816-11eb-17d8-b5509f40beca
md"""
### Common Denominator

**Common denominators**

You will have a list of rationals in the form:
 [ ($numer_1$, $denom_1$) , ... ($numer_n$, $denom_n$) ]

where all numbers are positive integers.

You have to produce a result in the form
 ($n_1$, $d$) ... ($n_n$, $d$), in which $d$ is as small as possible and
         $\frac{n_1}{d} \equiv \frac{numer_1}{denom_1}$ ... $\frac{n_n}{d} \equiv \frac{numer_n}{denom_n}$.

Example:

`convert_fracs([(1, 2), (1, 3), (1, 4)]) should be [(6, 12), (4, 12), (3, 12)]`

Note:

- Due to the fact that first translations were written long ago - more than 4 years - these translations have only irreducible fractions.
- Newer translations have some reducible fractions. To be on the safe side it is better to do a bit more work by simplifying fractions even if they don't have to be.
"""

# ╔═╡ 69309f32-6817-11eb-0e36-2788da6513af
const VT{T} = Vector{Tuple{T, T}} where T

# ╔═╡ 7bdcddec-6817-11eb-09a5-53e299b8771f
begin

function convert_fracs(frac_lst::VT{T})::VT{T} where T <: Integer
    @assert size(frac_lst, 1) > 0

    ## 1 - simplify all fractions if possible
    frac_lst = map(t -> _simplify(t...), frac_lst)

    size(frac_lst, 1) == 1 && return frac_lst

    ## 2 - determine common denominator
    den_prod = find_common_den(frac_lst)

    ## 3 - build new fractions
	[(num * (den_prod / den), den_prod) for (num, den) in frac_lst]
end

function _simplify(num::T, den::T)::Tuple{T, T} where T <: Integer
	gcd = _gcd(num, den)
	(num ÷ gcd, den ÷ gcd)
end

"""
greatest common divisor

ex.
    a=24, b=36 => 12
    a=2, b=8 => 2
    a=3, b=2 => 1
    a=9, b=6 => 3
"""
function _gcd(a::T, b::T)::T where T <: Integer
	a, b =  a > b ? (a, b) : (b, a)
	while true
        r = a % b
		r == 0 && (return b)
		r == 1 && (return r)
		a, b = (b, r)
    end
end

"""
lowest common multiple

ex.
    a=4, b=3 => 12
    a=2, b=8 => 8
    a=9 (3*3), b=6 (2*3) => 18
"""
function _lcm(a::T, b::T)::T where T <: Integer
    (a * b) ÷ _gcd(a, b)
end

function find_common_den(frac_lst::VT{T})::T where T <: Integer
    den_ary = [den for (_, den) in frac_lst]
    com_den = _lcm(den_ary[1], den_ary[2])

    if size(frac_lst, 1) > 2
		com_den = foldl((com_den, cur_den) -> _lcm(com_den, cur_den),
			den_ary[3:end];
			init=com_den
		)
    end
    return com_den
end

end

# ╔═╡ 07856cdc-681b-11eb-28d3-3b7d32a3afd1
begin
	@test _gcd(3, 4) == 1
	@test _gcd(4, 3) == 1
	@test _gcd(24, 36) == 12
	@test _gcd(8, 2) == 2
	@test _gcd(3, 2) == 1
	@test _gcd(9, 6) == 3
end

# ╔═╡ 1e4f2bce-681b-11eb-00cb-1baae1e5d600
begin
	@test _lcm(3, 4) == 12
	@test _lcm(4, 3) == 12

 	@test _lcm(1, 6) == 6
 	@test _lcm(6, 1) == 6

	@test _lcm(6, 3) == 6
	@test _lcm(3, 5) == 15
	@test _lcm(5, 3) == 15

 	@test _lcm(2, 7) == 14
 	@test _lcm(7, 2) == 14

 	@test _lcm(2*3, 2*2) == 12
end

# ╔═╡ 4c877dc8-681b-11eb-2f8d-bb469d56f840
with_terminal() do
	println(convert_fracs([(1, 2), (1, 3), (1, 4)]))
end

# ╔═╡ 952b7954-681c-11eb-0d11-b7ba99db8b64
begin
	@test convert_fracs([(2, 4)]) == [(1, 2)]
	@test convert_fracs([(1, 2), (1, 3), (1, 4)]) == [(6, 12), (4, 12), (3, 12)]
	@test convert_fracs([(2, 4), (2, 6), (2, 8)]) == [(6, 12), (4, 12), (3, 12)]

	@test convert_fracs([(1, 2), (1, 3)]) == [(3, 6), (2, 6)]

	@test convert_fracs([(1, 3), (1, 2)]) == [(2, 6), (3, 6)]
end

# ╔═╡ Cell order:
# ╟─d24d5ea0-6816-11eb-17d8-b5509f40beca
# ╠═6764c13e-681b-11eb-3527-77caf9d74f5b
# ╠═63a5111a-6817-11eb-2151-fb825d7879ba
# ╠═69309f32-6817-11eb-0e36-2788da6513af
# ╠═7bdcddec-6817-11eb-09a5-53e299b8771f
# ╠═07856cdc-681b-11eb-28d3-3b7d32a3afd1
# ╠═1e4f2bce-681b-11eb-00cb-1baae1e5d600
# ╠═4c877dc8-681b-11eb-2f8d-bb469d56f840
# ╠═952b7954-681c-11eb-0d11-b7ba99db8b64

max_stars_repo_path: Julia Kata/05_common_denominator.jl

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_path: Julia Kata/05_common_denominator.jl

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_path: Julia Kata/05_common_denominator.jl

max_forks_repo_name: pascal-p/julia-notebooks


text:
export Day2, Day2_1, Day2_2
abstract type Day2 <: Day end
abstract type Day2_1 <: Day2 end
abstract type Day2_2 <: Day2 end

struct Day2Line
    low::Int
    high::Int
    char::Char
    password::String
end

line_parser(::Type{<:Day2}, line) = begin
    (first, middle, last) = line |> strip |> split
    low, high = first |> s -> split(s, "-") .|> s -> parse(Int,s)
    Day2Line(low, high, middle[1], last)
end

solve(::Type{T}, lines::Vector) where T <: Day2 = lines |> lines -> valids(T, lines) |> sum

valids(::Type{T}, lines::Array{Day2Line,1}) where T <: Day2 = lines .|> l -> validpassword(T, l)

function validpassword(::Type{Day2_1}, line::Day2Line)
    count = 0
    for c in line.password
        if c == line.char
            count += 1
        end
    end
    line.low <= count <= line.high
end

function validpassword(::Type{Day2_2}, line::Day2Line)
    chars = Set([line.password[line.low], line.password[line.high]])
    (line.char in chars) & (length(chars) == 2)
end


text:
#ENV["CPLEX_STUDIO_BINARIES"]="/home/software/cplex/128/cplex/bin/x86-64_linux"

######################################################################################################################################################
#### Note: you will need to manually change the miop_formulation file to set a time limit of 120s before running this code, to reproduce output.
######################################################################################################################################################

using DelimitedFiles, LinearAlgebra, Random
include("src/core.jl")
function getProblemDataBeasley_minreturn(pathName)
    testData=readdlm(pathName)
    n=testData[1,1]
    μ=testData[2:n+1, 1]
    μ=convert(Array{Float64}, μ)
    σ=testData[2:n+1, 2]
    Σ=zeros(n,n)
    for iter=n+2:size(testData,1)
        # This loop inserts diagonal elements twice. This is ok as benign, and a "fix" would require more lines of code.
        Σ[testData[iter, 1], testData[iter, 2]]=σ[testData[iter, 1]]*σ[testData[iter, 2]]*testData[iter, 3]
        Σ[testData[iter, 2], testData[iter, 1]]=σ[testData[iter, 1]]*σ[testData[iter, 2]]*testData[iter, 3]
    end
    X=sqrt(Σ)
    k=5
    γ=100/sqrt(n)*ones(n)

    # Compute r_min and r_max, in order to add minimum return constraint.

    # r_min
    m=Model(solver=MosekSolver(MSK_IPAR_LOG=0, MSK_IPAR_MAX_NUM_WARNINGS=0));
    @variable(m, x[1:n]>=0.0);
    @constraint(m, sum(x)==1.0);
    @objective(m, Min,  ((1.0/(2.0*γ[1]))*x'*x +0.5*x'*Σ*x)');
    solve(m)
    r_min=μ'*getvalue(x)
    # r_max
    m2=Model(solver=MosekSolver(MSK_IPAR_LOG=0, MSK_IPAR_MAX_NUM_WARNINGS=0))
    @variable(m2, x2[1:n]>=0.0)
    @constraint(m2, sum(x2)==1.0)
    @objective(m2, Max,  (μ'*x2-(1.0/(2.0*γ[1]))*x2'*x2)')
    solve(m2)
    r_max=μ'*getvalue(x2)
    r_bar=r_min+0.3*(r_max-r_min)

    A=μ'
    l=[r_bar]
    u=[1e12]
    min_investment=zeros(n)
    f=size(X, 2)
    m=size(A, 1)
    # Constraining minimum return instead, so these entries all become 0.
    μ=zeros(n)
    Y=(X*X')\(X*μ)
    d=X'*Y-μ
    d=round.(d, digits=8) #avoids round-off errors when Σ is full rank

    return SparsePortfolioData(μ, Y, d, X, A, l, u, k, n, m, f, min_investment, γ)
end

for portfolioNumber=2:4
     pathName = "/Users/Ryancw2/Dropbox (MIT)/Sparse Portfolios Project/OptimalMarkowitzPortfolios.jl/data/port"
     pathName=pathName*string(portfolioNumber)*".txt"
     thePortfolioData=getProblemDataBeasley_minreturn(pathName)

     sparse_Portfolio=cutting_planes_portfolios(thePortfolioData; useSOCPLB=false, minReturnConstraint=true, useHeuristic=false);
     @show @test(abs.(sum(sparse_Portfolio.w)-1.0)<1e-4)
     println("Variance is:", sparse_Portfolio.portfolioVariance)
     println("Return is:", sparse_Portfolio.expectedReturn)
     println("Reg term is:", 0.5*sum(sparse_Portfolio.w[i]^2/thePortfolioData.γ[sparse_Portfolio.indices[i]] for i=1:size(sparse_Portfolio.indices,1)))
     open("results/Table5ResultsOutput.csv","a") do fp
     println(fp, "Instance, cardinality, regularizer, method, solvetime, gap, nodecount, iters")
     println(fp, "beasley"*string(portfolioNumber), ",", thePortfolioData.k, ",",
                       "100/sqrt(n)", ",", "OA", ",",  sparse_Portfolio.solveTime, ",", sparse_Portfolio.boundGap, ",", sparse_Portfolio.nodeCount, ",", sparse_Portfolio.numIters)
     end
     # You will need to manually check there are 50 (not 200) cuts applied here.
     sparse_Portfolio=cutting_planes_portfolios(thePortfolioData; useSOCPLB=false, minReturnConstraint=true, useHeuristic=false, usingKelleyPrimal=true, maxCutCallbacks=0);
     @show @test(abs.(sum(sparse_Portfolio.w)-1.0)<1e-4)
     println("Variance is:", sparse_Portfolio.portfolioVariance)
     println("Return is:", sparse_Portfolio.expectedReturn)
     println("Reg term is:", 0.5*sum(sparse_Portfolio.w[i]^2/thePortfolioData.γ[sparse_Portfolio.indices[i]] for i=1:size(sparse_Portfolio.indices,1)))
     open("results/Table5ResultsOutput_b.csv","a") do fp
     println(fp, "Instance, cardinality, regularizer, method, solvetime, gap, nodecount, iters")
     println(fp, "beasley"*string(portfolioNumber), ",", thePortfolioData.k, ",",
                       "100/sqrt(n)", ",", "OA-Kelley (50 cuts)", ",",  sparse_Portfolio.solveTime, ",", sparse_Portfolio.boundGap, ",", sparse_Portfolio.nodeCount, ",", sparse_Portfolio.numIters)
     end

     #Test against Big-M for correctness
     sparse_Portfolio_Cplex_bigM=cplex_raw_bigM(thePortfolioData, 120.0, 1e-8, 1)
     @show @test(abs.(sum(sparse_Portfolio_Cplex_bigM.w)-1.0)<1e-4)
     println("Variance is:", sparse_Portfolio_Cplex_bigM.portfolioVariance)
     println("Return is:", sparse_Portfolio_Cplex_bigM.expectedReturn)
     println("Reg term is:", 0.5*sum(sparse_Portfolio_Cplex_bigM.w[i]^2/thePortfolioData.γ[sparse_Portfolio_Cplex_bigM.indices[i]] for i=1:size(sparse_Portfolio_Cplex_bigM.indices,1)))
     open("results/Table5ResultsOutput.csv","a") do fp
     println(fp, "Instance, cardinality, regularizer, method, solvetime, gap, nodecount")
     println(fp, "beasley"*string(portfolioNumber), ",", thePortfolioData.k, ",",
                     "100/sqrt(n)", ",", "CPLEX Big-M", ",",  sparse_Portfolio_Cplex_bigM.solveTime, ",", sparse_Portfolio_Cplex_bigM.boundGap, ",", sparse_Portfolio_Cplex_bigM.nodeCount)
    end

    # #  # Test against MISOCP for correctness
    sparse_Portfolio_Cplex=cplex_raw_MISOCP(thePortfolioData, 120.0, 1e-8, 1)
    @show @test(abs.(sum(sparse_Portfolio_Cplex.w)-1.0)<1e-4)
    println("Variance is:", sparse_Portfolio_Cplex.portfolioVariance)
    println("Return is:", sparse_Portfolio_Cplex.expectedReturn)
    println("Reg term is:", 0.5*sum(sparse_Portfolio_Cplex.w[i]^2/thePortfolioData.γ[sparse_Portfolio_Cplex.indices[i]] for i=1:size(sparse_Portfolio_Cplex.indices,1)))
    open("results/Table5ResultsOutput.csv","a") do fp
    println(fp, "Instance, cardinality, regularizer, method, solvetime, gap, nodecount")
    println(fp, "beasley"*string(portfolioNumber), ",", thePortfolioData.k, ",",
                        "100/sqrt(n)", ",", "CPLEX MISOCP", ",",  sparse_Portfolio_Cplex.solveTime, ",", sparse_Portfolio_Cplex.boundGap, ",", sparse_Portfolio_Cplex.nodeCount)
    end
end


text:
# Binomial distribution

export Binomial
import Base
using SpecialFunctions

probit(p) = sqrt2 * erfinv(2p - 1)
Φ(z) = (1 + erf(invsqrt2 * z))/2

@parameterized Binomial(n,p) ≪ CountingMeasure(ℤ[0:∞])

(d::Binomial ≪ ::CountingMeasure{IntegerRange{a,b}}) where {a,b} = a ≤ 0 && b ≥ d.n

(::CountingMeasure{IntegerRange{a,b}} ≪ ::Binomial) where {a,b} = a ≥ 0 && b ≤ d.n

###############################################################################
@kwstruct Binomial(n, p)

function logdensity(d::Binomial{(:n, :p)}, y)
    (n, p) = (d.n, d.p)
    return -log1p(n) - logbeta(n - y + 1, y + 1) + xlogy(y, p) + xlog1py(n - y, -p)
end

function Base.rand(rng::AbstractRNG, T::Type, d::Binomial{(:n,:p)})
    rand(rng, Dists.Binomial(d.n, d.p))
end

###############################################################################
@kwstruct Binomial(n, logitp)

function logdensity(d::Binomial{(:n, :logitp)}, y)
    n = d.n
    x = d.logitp
    return  -log1p(n) - logbeta(n - y + 1, y + 1) + y * x - n * log1pexp(x)
end

function Base.rand(rng::AbstractRNG, d::Binomial{(:n,:logitp)})
    rand(rng, Dists.Binomial(d.n, logistic(d.logitp)))
end

###############################################################################
@kwstruct Binomial(n, probitp)

function logdensity(d::Binomial{(:n, :probitp)}, y)
    n = d.n
    z = d.probitp
    return  -log1p(n) - logbeta(n - y + 1, y + 1)  + xlogy(y, Φ(z)) + xlogy(n-y, Φ(-z))
end

function Base.rand(rng::AbstractRNG, d::Binomial{(:n, :probitp)})
    rand(rng, Dists.Binomial(d.n, Φ(d.probitp)))
end

distproxy(d::Binomial{(:n, :p)}) = Dists.Binomial(d.n, d.p)
distproxy(d::Binomial{(:n,:logitp)}) = Dists.Binomial(d.n, logistic(d.logitp))
distproxy(d::Binomial{(:n,:probitp)}) = Dists.Binomial(d.n, Φ(d.probitp))

asparams(::Type{<:Binomial}, ::Val{:p}) = as𝕀
asparams(::Type{<:Binomial}, ::Val{:logitp}) = asℝ
asparams(::Type{<:Binomial}, ::Val{:probitp}) = asℝ


text:
export svdl

import Base: size, getindex, Matrix
import LinearAlgebra.svd
using LinearAlgebra

"""
Matrix of the form
```
 d1          a_1
    d2       a_2
      ...    ...
         d_l a_l
             d_l+1 e_l+1
                   ...  e_k-1
                        d_k
```
"""
mutable struct BrokenArrowBidiagonal{T} <: AbstractMatrix{T}
    dv::Vector{T}
    av::Vector{T}
    ev::Vector{T}
end

size(B::BrokenArrowBidiagonal) = (n=length(B.dv); (n, n))

function size(B::BrokenArrowBidiagonal, n::Int)
    if n==1 || n==2
        return length(B.dv)
    else
        throw(ArgumentError("invalid dimension $n"))
    end
end

function getindex(B::BrokenArrowBidiagonal{T}, i::Int, j::Int) where {T}
    n = size(B, 1)
    k = length(B.av)
    if !(1 ≤ i ≤ n && 1 ≤ j ≤ n)
        throw(BoundsError())
    end

    if i == j
        return B.dv[i]
    elseif i ≤ k && j == k+1
        return B.av[i]
    elseif i > k && j == i+1
        return B.ev[i-k]
    else
        return zero(T)
    end
end

function Matrix(B::BrokenArrowBidiagonal{T}) where {T}
    n = size(B, 1)
    k = length(B.av)
    M = zeros(T, n, n)
    for i=1:n
        M[i, i] = B.dv[i]
    end
    for i=1:k
        M[i,k+1] = B.av[i]
    end
    for i=k+1:n-1
        M[i, i+1] = B.ev[i-k]
    end
    M
end

svd(B::BrokenArrowBidiagonal) = svd(Matrix(B)) #XXX This can be much faster

"""
Partial factorization object which is an approximation of a matrix

    A ≈ P * [B 0; 0 β] * Q
"""
mutable struct PartialFactorization{T, Tr} <: Factorization{T}
    P :: Matrix{T}
    Q :: Matrix{T}
    B :: AbstractMatrix{Tr}
    β :: Tr
end

####################
# API method calls #
####################

"""
    svdl(A) -> Σ, L, [history]

Compute some singular values (and optionally vectors) using Golub-Kahan-Lanczos
bidiagonalization [^Golub1965] with thick restarting [^Wu2000].

If `log` is set to `true` is given, method will output a tuple `X, L, ch`. Where
`ch` is a `ConvergenceHistory` object. Otherwise it will only return `X, L`.

# Arguments

- `A` : The matrix or matrix-like object whose singular values are desired.

## Keywords

- `nsv::Int = 6`: number of singular values requested;
- `v0 = random unit vector`: starting guess vector in the domain of `A`.
  The length of `q` should be the number of columns in `A`;
- `k::Int = 2nsv`: maximum number of Lanczos vectors to compute before restarting;
- `j::Int = nsv`: number of vectors to keep at the end of the restart.
  We don't recommend j < nsv;
- `maxiter::Int = minimum(size(A))`: maximum number of iterations to run;
- `verbose::Bool = false`: print information at each iteration;
- `tol::Real = √eps()`: maximum absolute error in each desired singular value;
- `reltol::Real=√eps()`: maximum error in each desired singular value relative to the
  estimated norm of the input matrix;
- `method::Symbol=:ritz`: restarting algorithm to use. Valid choices are:
  1. `:ritz`: Thick restart with Ritz values [Wu2000].
  2. `:harmonic`: Restart with harmonic Ritz values [Baglama2005].
- `vecs::Symbol = :none`: singular vectors to return.
  1. `:both`: Both left and right singular vectors are returned.
  2. `:left`: Only the left singular vectors are returned.
  3. `:right`: Only the right singular vectors are returned.
  4. `:none`: No singular vectors are returned.
- `dolock::Bool=false`: If `true`, locks converged Ritz values, removing them
  from the Krylov subspace being searched in the next macroiteration;
- `log::Bool = false`: output an extra element of type `ConvergenceHistory`
  containing extra information of the method execution.

# Return values

**if `log` is `false`**

- `Σ`: list of the desired singular values if `vecs == :none` (the default), otherwise
  returns an `SVD` object with the desired singular vectors filled in;
- `L`: computed partial factorizations of A.

**if `log` is `true`**

- `Σ`: list of the desired singular values if `vecs == :none` (the default),
otherwise returns an `SVD` object with the desired singular vectors filled in;
- `L`: computed partial factorizations of A;
- `history`: convergence history.

**ConvergenceHistory keys**

- `:betas` => `betas`: The history of the computed betas.
- `:Bs` => `Bs`: The history of the computed projected matrices.
- `:ritz` => `ritzvalhist`: Ritz values computed at each iteration.
- `:conv` => `convhist`: Convergence data.

[^Golub1965]:
    Golub, Gene, and William Kahan. "Calculating the singular values and pseudo-inverse
    of a matrix." Journal of the Society for Industrial and Applied Mathematics,
    Series B: Numerical Analysis 2.2 (1965): 205-224.
[^Wu2000]:
    Wu, Kesheng, and Horst Simon. "Thick-restart Lanczos method for large symmetric
    eigenvalue problems." SIAM Journal on Matrix Analysis and Applications 22.2
    (2000): 602-616.
"""
function svdl(A;
    nsv::Int=6, k::Int=2nsv, tol::Real=√eps(),
    maxiter::Int=minimum(size(A)), method::Symbol=:ritz, log::Bool=false, kwargs...
    )
    history = ConvergenceHistory(partial=!log)
    history[:tol] = tol
    reserve!(BitArray, history,:conv, maxiter)
    reserve!(history,[:ritz,:resnorm], maxiter, k)
    Bs_type = (method == :ritz) ? BrokenArrowBidiagonal : UpperTriangular
    reserve!(Bs_type, history,:Bs, maxiter)
    reserve!(history,:betas, maxiter)
    X, L = svdl_method!(history, A, nsv; k=k, tol=tol, maxiter=maxiter, method=method, kwargs...)
    log && shrink!(history)
    log ? (X, L, history) : (X, L)
end

#########################
# Method Implementation #
#########################

function svdl_method!(log::ConvergenceHistory, A, l::Int=min(6, size(A,1)); k::Int=2l,
    j::Int=l, v0::AbstractVector = Vector{eltype(A)}(randn(size(A, 2))) |> x->rmul!(x, inv(norm(x))),
    maxiter::Int=minimum(size(A)), tol::Real=√eps(), reltol::Real=√eps(),
    verbose::Bool=false, method::Symbol=:ritz, vecs=:none, dolock::Bool=false)

    T0 = time_ns()
    @assert k>1
    L = build(log, A, v0, k)

    iter = 0
    local F
    for iter in 1:maxiter
        nextiter!(log)

        #@assert size(L.B) == (k, k)
        F = svd(L.B)::LinearAlgebra.SVD{eltype(v0), typeof(real(one(eltype(v0)))), Matrix{eltype(v0)}}
        iter==1 && @assert eltype(F)==eltype(v0)
        if method == :ritz
            thickrestart!(A, L, F, j)
        elseif method == :harmonic
            harmonicrestart!(A, L, F, j)
        else
            throw(ArgumentError("Unknown restart method $method"))
        end
        extend!(log, A, L, k)
        if verbose
            elapsedtime = round((time_ns()-T0)*1e-9, digits=3)
            info("Iteration $iter: $elapsedtime seconds")
        end

        conv = isconverged(L, F, l, tol, reltol, log, verbose)

        push!(log, :conv, conv)
        push!(log, :ritz, F.S[1:k])
        push!(log, :Bs, deepcopy(L.B))
        push!(log, :betas, L.β)

        #Lock
        if method == :ritz && dolock
            for i in eachindex(conv)
                if conv[i]
                    L.B.av[i] = 0
                end
            end
        end
        all(conv) && (setconv(log, true); break)
    end
    shrink!(log)

    #Compute singular vectors as necessary and return them in the output
    values = F.S[1:l]
    m, n = size(A)

    leftvecs = if vecs == :left || vecs == :both
        L.P*view(F.U, :, 1:l)
    else
        zeros(eltype(v0), m, 0)
    end

    rightvecs = if vecs == :right || vecs == :both
        (view(L.Q, :, 1:size(L.Q,2)-1)*view(F.V, :, 1:l))'
    else
        zeros(eltype(v0), 0, n)
    end

    if vecs == :none
        values, L
    else
        LinearAlgebra.SVD(leftvecs, values, rightvecs), L
    end
end

"""
    isconverged(L, F, k, tol, reltol, log, verbose=false)

Determine if any singular values in a partial factorization have converged.

# Arguments

- `L::PartialFactorization`: a `PartialFactorization` computed by an iterative
method such as `svdl`;
- `F::LinearAlgebra.SVD`: a `SVD` factorization computed for `L.B`;
- `k::Int` : number of singular values to check;
- `tol::Real`: absolute tolerance for a Ritz value to be considered converged;
- `reltol::Real`: relative tolerance for a Ritz value to be considered converged;
- `verbose::Bool = false`: if `true`, prints out all the results of convergence tests.

# Implementation note

This convergence test routine uses a variety of different tests.
1. The crudest estimate of the error bound is a simple error bound which dates
back at least to Wilkinson's classic book [^Wilkinson1965] Ch.3 §53 p.170.
2. When the Ritz values become sufficiently well-separated, more refined estimates
can be derived from the Rayleigh-Ritz properties of the Krylov process, as
described in [^Wilkinson1965] Ch.3 §54-55 p.173, [^Yamamoto1980], [^Ortega1990],
[^Geurts1982] and [^Deif1989].

[^Wilkinson1965]:
    Wilkinson, James Hardy. The algebraic eigenvalue problem.
    Vol. 87. Oxford: Clarendon Press, 1965.
[^Yamamoto1980]:
    Yamamoto, Tetsuro. "Error bounds for computed eigenvalues
    and eigenvectors." Numerische Mathematik 34.2 (1980): 189-199.
[^Ortega1990]:
    Ortega, James M. Numerical analysis: a second course.
    Society for Industrial and Applied Mathematics, 1990.
[^Geurts1982]:
    Geurts, A. J. "A contribution to the theory of condition."
    Numerische Mathematik 39.1 (1982): 85-96.
[^Deif1989]:
    Deif, A. "A relative backward perturbation theorem for the eigenvalue
    problem." Numerische Mathematik 56.6 (1989): 625-626.
"""
function isconverged(L::PartialFactorization, F::LinearAlgebra.SVD, k::Int, tol::Real,
    reltol::Real, log::ConvergenceHistory, verbose::Bool=false
    )

    @assert tol ≥ 0

    σ = F.S[1:k]
    Δσ= L.β * abs.(F.U[end, 1 : k])

    #Best available eigenvalue bounds
    δσ = copy(Δσ)

    #Update better error bounds from Rayleigh-Ritz
    #In theory these only apply if we know all the values
    #But so long as the gap between the converged Ritz values and all the
    #others is larger than the smallest empirical spectral gap d then we're fine.
    #Can also get some eigenvector statistics!!
    if k > 1
        d = Inf
        for i in eachindex(σ), j=1:i-1
            d = min(d, abs(σ[i] - σ[j]))
        end
        verbose && println("Smallest empirical spectral gap: ", d)
        #Chatelein 1993 - normwise backward error associated with approximate invariant subspace
        #Use the largest singular (Ritz) value to estimate the 2-norm of the matrix
        verbose && println("Normwise backward error associated with subspace: ", L.β/σ[1])
        for i in eachindex(Δσ)
            α = Δσ[i]

            #Simple error bound
            if 2α ≤ d
                #Rayleigh-Ritz bounds
                x = abs(α/(d-α)*√(1+(α/(d-α))^2))
                verbose && println("Rayleigh-Ritz error bound on eigenvector: $x")
                #δσ[i] = min(δσ[i], x)

                y = α^2/d #[Wilkinson:Ch.3 Appendix (4), p.188]
                verbose && println("Rayleigh-Ritz error bound on eigenvalue: $y")
                δσ[i] = min(δσ[i], y)
            end

	    verbose && println("Ritz value ", i, ": ", σ[i], " ± ", round(δσ[i], sigdigits=3))

            #Estimate of the normwise backward error [Deif 1989]
            ##Use the largest singular (Ritz) value to estimate the 2-norm of the matrix
            verbose && println("Normwise backward error estimate: ", α/σ[1])

            #Geurts, 1982 - componentwise backward error also known
            #A. J. Geurts, (1982), A contribution to the theory of condition,
            #Numer.  Math., #39, 85-96.
        end
    end

    #Estimate condition number and see if two-sided reorthog is needed
    if verbose && (F.S[1]/F.S[end]) > 1/√eps()
        @warn("Two-sided reorthogonalization should be used but is not implemented")
    end

    push!(log, :resnorm, δσ[1:k])
    conv = (δσ[1:k] .< max(tol, reltol*σ[1]))::BitVector
end

#Hernandez2008
function build(log::ConvergenceHistory, A, q::AbstractVector{T}, k::Int) where {T}
    m, n = size(A)
    Tr = typeof(real(one(T)))
    β = norm(q)
    q .*= inv(β)
    p = A*q
    α = convert(Tr, norm(p))
    p .*= inv(α)
    bidiag = Bidiagonal([α], Tr[], :U)
    extend!(log, A, PartialFactorization(reshape(p, m, 1), reshape(q, n, 1), bidiag, β), k)
end

"""
    thickrestart!(A, L,F, l)

Thick restart (with ordinary Ritz values)

# References

[^Hernandez2008]

"""
function thickrestart!(A, L::PartialFactorization{T,Tr},
        F::LinearAlgebra.SVD{Tr,Tr}, l::Int) where {T,Tr}

    k = size(F.V, 1)
    m, n = size(A)
    #@assert size(L.P) == (m, k)
    #@assert size(L.Q) == (n, k+1)

    Q = view(L.Q, :,1:k)*view(F.V, :,1:l)
    L.Q = [Q view(L.Q, :, k+1)]
    #Be pedantic about ensuring normalization
    #L.Q = qr(L.Q)[1]
    #@assert all([norm(L.Q[:,i]) ≈ 1 for i=1:size(L.Q,2)])

    f = A*view(L.Q, :, l+1)
    ρ = L.β * reshape(F.U[end, 1:l], l)
    L.P = view(L.P, :, 1:k)*view(F.U, :, 1:l)

    #@assert ρ[i] ≈ f⋅L.P[:, i]
    f -= L.P*ρ
    α = convert(Tr, norm(f))
    f .*= inv(α)
    L.P = [L.P f]

    g = A'f - α*L.Q[:, end]
    L.β = β = norm(g)
    L.B = BrokenArrowBidiagonal([F.S[1:l]; α], ρ, typeof(β)[])
    #@assert size(L.P, 2) == size(L.Q, 2) == size(L.B, 2)
    L
end

"""
    harmonicrestart!(A, L, F, k) -> L

Thick restart with harmonic Ritz values. See [^Baglama2005] but note that
they have P and Q swapped relative to our notation, which follows that
of [^Hernandez2008]

[^Baglama2005]:
    Baglama, James, and Lothar Reichel. "Augmented implicitly restarted
    Lanczos bidiagonalization methods." SIAM Journal on Scientific
    Computing 27.1 (2005): 19-42.
[^Hernandez2008]:
    Vicente Hernández, José E. Román, and Andrés Tomás. "A robust and
    efficient parallel SVD solver based on restarted Lanczos bidiagonalization."
    Electronic Transactions on Numerical Analysis 31 (2008): 68-85.

"""
function harmonicrestart!(A, L::PartialFactorization{T,Tr},
        F::LinearAlgebra.SVD{Tr,Tr}, k::Int) where {T,Tr}

    m = size(L.B, 1)::Int
    #@assert size(L.P,2)==m==size(L.Q,2)-1

    F0 = F# svd(L.B)
    ρ = L.β*F0.U[end,:] #Residuals of singular values

    #Construct broken arrow matrix
    BA = [Matrix(Diagonal(F0.S)) ρ]
    F2 = svd!(BA; full=true)

    #Take k largest triplets
    Σ = (F2.S::Vector{Tr})[1:k]
    U = F0.U*view(F2.U,:,1:k)
    M = Matrix{T}(I, m+1, m+1)
    M[1:m, 1:m] = F0.V::Union{Matrix{T}, Adjoint{T,Matrix{T}}}
    M = M * F2.V
    Mend = M[end, 1:k]
    #Compute scaled residual from the harmonic Ritz problem
    r0 = zeros(Tr, m)
    r0[end] = 1
    if isa(L.B, Bidiagonal)
        @assert L.B.uplo == 'U'
    end
    r = try
        #(L.B\r0)
        ldiv!(L.B, r0)
    catch exc
        if isa(exc, LinearAlgebra.LAPACKException) ||
            isa(exc, LinearAlgebra.SingularException) #B\r is singular
            pinv(Matrix(L.B))*r0
        else rethrow(exc) end
    end::Vector{Tr}
    r .*= L.β
    M::Matrix{T} = view(M,1:m, :) + r*view(M,m+1:m+1,:)

    M2 = zeros(T, m+1, k+1)
    M2[1:m, 1:k] = M[:,1:k]
    M2[1:m, k+1] = -r
    M2[m+1, k+1] = 1
    QRF = qr(M2)
    Q, R = QRF.Q, QRF.R

    Q = L.Q*view(Q,:,1:k+1)
    P = L.P*view(U,:,1:k)

    R = view(R,1:k+1,1:k) + view(R,:,k+1)*Mend'

    f = A*view(Q,:,k+1)
    f -= P*(P'f)
    α = convert(Tr, norm(f))
    f .*= inv(α)
    P = [P f]
    B = UpperTriangular{Tr,Matrix{Tr}}([(Diagonal(Σ)*triu(R')); zeros(Tr,1,k) α])
    #@assert size(P, 2) == size(B, 1) == size(Q, 2)
    g = A'f
    q = view(Q,:,k+1)
    #g-= (g⋅q)*q
    g .-= (g⋅q)*q
    β = convert(Tr, norm(g))
    g .*= inv(β)
    #@assert size(P, 2) == size(Q, 2) == size(B, 2)
    L.P = P
    L.Q = Q
    L.B = B
    L.β = β
    L
end

"""
    extend!{T,Tr}(A, L, k orthleft, orthright, α)

Extend a PartialFactorization L using GKL bidiagonalization with k extra pairs
of Lanczos vectors.

# Arguments

- `A`: matrix or linear map generating the Lanczos vectors;
- `L::PartialFactorization`: partial factorization;
- `orthleft::Bool = false`: orthogonalize left Lanczos vectors;
- `orthright::Bool = true`: orthogonalize right Lanczos vectors;
- `α::Real = 1/√2`: criterion for doing a second reorthogonalization.

## Implementation details

The implementation mostly follows the description in [^Simon2000] and [^Hernandez2008].

The reorthogonalization method used is using double classical Gram-Schmidt
full reorthogonalization. As explained in the numerical analysis literature by
Kahan, Golub, Rutishauser, and others in the 1970s, double classical
Gram-Schmidt reorthogonalization always suffices to keep vectors orthogonal to
within machine precision. As described in [^Bjorck2015], `α` is a threshold
for determinining when the second orthogonalization is necessary. -log10(α) is
the number of (decimal) digits lost due to cancellation. Common choices are
`α=0.1` [Rutishauser] and `α=1/√2` [Daniel1976] (our default).

In most situations it suffices to orthogonalize either the left vectors or the
right vectors, except when the matrix norm exceeds `1/√eps(eltype(A))`, in
which case it will be necessary to orthogonalize both sets of vectors. See
[^Simon2000].

[^Bjorck2015]:
    Björck, Åke. Numerical methods in matrix computations.
    New York: Springer, 2015.
[^Simon2000]:
    Simon, Horst D., and Hongyuan Zha. "Low-rank matrix approximation
    using the Lanczos bidiagonalization process with applications."
    SIAM Journal on Scientific Computing 21.6 (2000): 2257-2274.

[^Daniel1976]:
    Daniel, James W., et al. "Reorthogonalization and stable algorithms
    for updating the Gram-Schmidt QR factorization." Mathematics of
    Computation 30.136 (1976): 772-795.
```
"""
function extend!(
    log::ConvergenceHistory, A, L::PartialFactorization{T, Tr}, k::Int,
    orthleft::Bool=false, orthright::Bool=true, α::Real = 1/√2
    ) where {T,Tr}

    l = size(L.B, 2)::Int-1
    p = L.P[:,l+1]

    m, n = size(A)
    q = zeros(T, n)
    #@assert l+1 == size(L.B, 1) == size(L.B, 2) == size(L.P, 2)

    if isa(L.B, UpperTriangular) #Cannot be appended
        B = zeros(Tr, k, k)
        B[1:size(L.B,1), 1:size(L.B,2)] = L.B
        L.B = UpperTriangular(B)
        #@assert size(L.B) == (k, k)
    end

    β = L.β
    for j=l+1:k
        log.mtvps+=1
        mul!(q, adjoint(A), p) #q = A'p

        if orthright #Orthogonalize right Lanczos vector
            #Do double classical Gram-Schmidt reorthogonalization
            oldqnorm = norm(q)
            q -= L.Q*(L.Q'q)
            if norm(q) ≤ α * oldqnorm
                q -= L.Q*(L.Q'q)
            end
        end

        β = norm(q)
        q .*= inv(β)

        L.Q = [L.Q q]
        j==k && break

        log.mvps+=1
        #p = A*q - β*p
        mul!(p, A, q)
        p .-= β*view(L.P, :, j)

        if orthleft #Orthogonalize left Lanczos vector
            #Do double classical Gram-Schmidt reorthogonalization
            oldpnorm = norm(p)
            p -= L.P*(L.P'p)
            if norm(p) ≤ α * oldpnorm
                p -= L.P*(L.P'p)
            end
        end

        α = norm(p)
        p .*= inv(α)
        if isa(L.B, Bidiagonal) || isa(L.B, BrokenArrowBidiagonal)
            push!(L.B.dv, α)
            push!(L.B.ev, β)
        else
            L.B[j+1, j+1] = α
            L.B[j  , j+1] = β
        end
        L.P = [L.P p]
    end
    L.β = β
    L
end

let
    B = BrokenArrowBidiagonal([1, 2, 3], [1, 2], Int[])
    @assert Matrix(B) == [1 0 1; 0 2 2; 0 0 3]
end


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using CairoMakie, ESDL
CairoMakie.activate!() # HIDE
let
    cube = esdc()
    cubeTemP =subsetcube(cube,variable=["air_temperature_2m", "precipitation"],
        time = Date(2003,1,1))
    #data_preci = readcubedata(cubeTemP)
    #datpreci = data_preci.data
    lon = getAxis("lon", cubeTemP).values
    lat = getAxis("lat", cubeTemP).values
    #time = getAxis("time", cube).values;
    precipitation = cubeTemP[var="precipitation"].data[:,:]
    temperature = cubeTemP[var="air_temperature_2m"].data[:,:]

    fig = Figure(resolution = (1250, 1200), font = "sans", fontsize = 24)
    ax1 = Axis(fig, xlabel = "lon", ylabel = "lat")
    ax2 = Axis(fig, xlabel = "lon", ylabel = "lat")

    map1 = heatmap!(ax1, lon, lat, temperature, colormap = :plasma)
    cbar1 = Colorbar(fig, map1, label = "temperature", ticklabelsize = 18,
        labelpadding = 5, width = 10)
    map2 = heatmap!(ax2, lon, lat, precipitation, colormap = :turbo,
        colorrange = (5,30))
    cbar2 = Colorbar(fig, map2, label = "precipitation", ticklabelsize = 18,
            labelpadding = 5, width = 10) #scale = log10 not working properly

    ax1.xticks = -180:60:180
    ax1.yticks = -90:30:90
    ax1.xtickformat = "{:d}ᵒ"
    ax1.ytickformat = "{:d}ᵒ"
    ax2.xticks = -180:60:180
    ax2.yticks = -90:30:90
    ax2.xtickformat = "{:d}ᵒ"
    ax2.ytickformat = "{:d}ᵒ"
    hidespines!(ax1, :t, :r)
    hidexdecorations!(ax1, ticks = false)
    hidespines!(ax2, :t, :r)
    fig[1,1] = ax1
    fig[1,2] = cbar1
    fig[2,1] = ax2
    fig[2,2] = cbar2
    fig
    save(joinpath(@__DIR__, "output", "temPrecESDL.png"), fig, px_per_unit = 2.0) # HIDE
end

using Pkg # HIDE
Pkg.status(["CairoMakie","ESDL"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
# n = mult(3, 4) #=> ErrorException("mult not defined"))
function mult(x, y)
    println("x is $x and y is $y")
    x * y
end

n = mult(3, 4) #=> 12
# one-line assignment syntax:
mult(x, y) = x * y
f(x,y) = x^3-y+x*y
println(f(3,2) ) # 31

function mult(x, y)
    println("x is $x and y is $y")
    if x == 1 return y end
    x * y
end

function mult(x::Float64, y::Float64)
    println("x is $x and y is $y")
    x * y
end
# following gives error:
# mult(5, 6)
# ERROR: `mult` has no method matching mult(::Int64, ::Int64)

function multi(n,m)
    n*m, div(n,m), n%m
end
x, y, z = multi(8, 2)

function varargs(n, m, args...)
    println("arguments: $n $m $args")
end
varargs(1, 2, 3, 4) # n is 1, m is 2, and args has value (3, 4)

x = (3, 4)
varargs2(1, 2, x...) # args is (1,2,3,4)
x = [10, 11, 12]
varargs2(1, 2, x...) # args is (1,2,10,11,12)
varargs2(k1="name1", k2="name2, k3=7)

apply(mult, 3, 4) #=> 12

∑(x,y) = x + y
∑(3, 4) #=> 7

# anonymous functions:
function (x)
    x + 2
end

g = (x, y) -> x^3-y+x*y
println(g(3,2) ) # 31

x -> x + 2
() -> x * 5

3+4 #=> 7
+(3,4) #=> 7

m = mult
m(6, 6) #=> 36

max_stars_repo_path: Module 1/Chapter03/Old/functions101.jl

max_issues_repo_path: Module 1/Chapter03/Old/functions101.jl

max_forks_repo_path: Module 1/Chapter03/Old/functions101.jl


text:
# Command: Perform Actions
# """
#     actions!(session::Session, actions::AbstractVector)
# """
# function actions!(session::Session, actions::AbstractVector)
#     @unpack addr, id = session
#     x = try
#     response = HTTP.post("$addr/session/$id/actions",
#                          [("Content-Type" => "application/json")],
#                          JSON3.write("parameters" => "moveto!Element($(hoverOver.id))"))
#     catch err
#         err
#     end
#     @assert response.status == 200
#     nothing
# end


text:
using Pkg
Pkg.activate(normpath(joinpath(@__DIR__, "..")))

using DataFrames
using CSV
using ArgParse
using Logging
using Printf
using ATools
include("../src/ANema.jl")

import ANema: getlorpath

logger = SimpleLogger(stdout, Logging.Warn)
old_logger = global_logger(logger)

function makelors(args)

	rdir   = args["dir"]
	conf   = args["conf"]
	pathf = ANema.pathfinder(rdir, conf)
	println("pathfinder tuple =", pathf)

	n3p = DataFrame(CSV.File(pathf.phipath))
	n3z = DataFrame(CSV.File(pathf.zpath))

	n3mp = ANema.setunits(n3p)
	n3mz = ANema.setunits(n3z)

	println("phi data frame has length of ", nrow(n3mp))
	println("z data frame has length of ", nrow(n3mz))
	println("phi data frame names ", names(n3mp))
	println("z data frame names ", names(n3mz))

	lors_first_true = ANema.dftolor(n3mz, ANema.dtfirst, ANema.postrue)
	ANema.write_lors_hdf5(getlorpath(pathf.lordir,pathf.zlp, "first-true"), lors_first_true)
	lors_minimum_reco = ANema.dftolor(n3mz, ANema.dtminimum, ANema.posreco)
	ANema.write_lors_hdf5(getlorpath(pathf.lordir,pathf.zlp, "minimum-reco"), lors_minimum_reco)
	lors_minimum_recall = ANema.dftolor(n3mz, ANema.dtminimum, ANema.posrecall)
	ANema.write_lors_hdf5(getlorpath(pathf.lordir,pathf.zlp, "minimum-recall"), lors_minimum_recall)
	lors_average_recall = ANema.dftolor(n3mz, ANema.dtaverage, ANema.posrecall)
	ANema.write_lors_hdf5(getlorpath(pathf.lordir,pathf.zlp, "average-recall"), lors_average_recall)

	lors_first_true = ANema.dftolor(n3mp, ANema.dtfirst, ANema.postrue)
	ANema.write_lors_hdf5(getlorpath(pathf.lordir,pathf.philp, "first-true"), lors_first_true)
	lors_minimum_reco = ANema.dftolor(n3mp, ANema.dtminimum, ANema.posreco)
	ANema.write_lors_hdf5(getlorpath(pathf.lordir,pathf.philp, "minimum-reco"), lors_minimum_reco)
	lors_minimum_recall = ANema.dftolor(n3mp, ANema.dtminimum, ANema.posrecall)
	ANema.write_lors_hdf5(getlorpath(pathf.lordir,pathf.philp, "minimum-recall"), lors_minimum_recall)
	lors_average_recall = ANema.dftolor(n3mp, ANema.dtaverage, ANema.posrecall)
	ANema.write_lors_hdf5(getlorpath(pathf.lordir,pathf.philp, "average-recall"), lors_average_recall)

end

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--dir", "-d"
            help = "directory with nema data"
            arg_type = String
            default = "../data"
		"--conf", "-c"
            help = "configuration"
            arg_type = String
    end

    return parse_args(s)
end

function main()
	parsed_args = parse_commandline()
	println("Running makelors with arguments", parsed_args)
	makelors(parsed_args)
end

@time main()

max_stars_repo_name: jjgomezcadenas/mytest

max_issues_repo_name: jjgomezcadenas/mytest

max_forks_repo_name: jjgomezcadenas/mytest


text:
"""
    WeightedCovMatrix(T = Float64)

Weighted covariance matrix, tracked as a matrix of type `T`.

*After* a call to `cov` the covariance matrix is stored in `o.C`.

# Example:
    o = fit!(WeightedCovMatrix(), rand(100, 4), rand(100))
    sum(o)
    mean(o)
    var(o)
    std(o)
    cov(o)
    cor(o)
"""
mutable struct WeightedCovMatrix{T} <: WeightedOnlineStat{VectorOb}
    C::Matrix{T}
    A::Matrix{T}
    b::Vector{T}
    W::T
    W2::T
    n::Int
    function WeightedCovMatrix{T}(
            C = zeros(T, 0, 0),
            A = zeros(T, 0, 0),
            b = zeros(T, 0),
            W = T(0),
            W2 = T(0),
            n = 0
        ) where T
        new{T}(C, A, b, W, W2, n)
    end
end

function WeightedCovMatrix(
        C::Matrix{T},
        A::Matrix{T},
        b::Vector{T},
        W::T,
        W2::T,
        n::Int
    ) where T
    WeightedCovMatrix{T}(C, A, b, W, W2, n)
end

WeightedCovMatrix(::Type{T}, p::Int=0) where T =
    WeightedCovMatrix(zeros(T, p, p), zeros(T, p, p), zeros(T, p), T(0), T(0), 0)
WeightedCovMatrix() = WeightedCovMatrix(Float64)

function _fit!(o::WeightedCovMatrix{T}, x, w) where T
    xx = convert(Vector{T}, x)
    ww = convert(T, w)

    o.n += 1
    γ1 = T(1) / o.n
    o.W = smooth(o.W, ww, γ1)
    o.W2 = smooth(o.W2, ww*ww, γ1)
    γ2 = ww / (o.W * o.n)
    if isempty(o.A)
        p = length(xx)
        o.b = zeros(T, p)
        o.A = zeros(T, p, p)
        o.C = zeros(T, p, p)
    end
    smooth!(o.b, xx, γ2)
    smooth_syr!(o.A, xx, γ2)
end

function _fit!(o::WeightedCovMatrix{T1},
               x::AbstractVector{Union{T2, Missing}}, w) where {T1, T2}
    if !mapreduce(ismissing, |, x)
        xx = convert(Vector{T1}, x)
        _fit!(o, xx, w)
    end
    return o
end
_fit!(o::WeightedCovMatrix, x, w::Missing) = o

function _merge!(o::WeightedCovMatrix{T}, o2::WeightedCovMatrix) where T
    o2_A = convert(Matrix{T}, o2.A)
    o2_b = convert(Vector{T}, o2.b)
    o2_W = convert(T, o2.W)
    o2_W2 = convert(T, o2.W2)

    if isempty(o.A)
        o.C = convert(Matrix{T}, o2.C)
        o.A = o2_A
        o.b = o2_b
        o.W = o2_W
        o.W2 = o2_W2
        o.n = o2.n
    else
        n = o.n + o2.n
        W = smooth(o.W, o2_W, o2.n / n)
        γ = (o2_W * o2.n) / (W * n)
        smooth!(o.A, o2_A, γ)
        smooth!(o.b, o2_b, γ)

        o.n = n
        o.W = W
        o.W2 = smooth(o.W2, o2_W2, o2.n / o.n)
    end

    return o
end

nvars(o::WeightedCovMatrix) = size(o.A, 1)

function OnlineStatsBase.value(o::WeightedCovMatrix)
    # o.A is only the upper triangle:
    # o.C .= o.A .- o.b .* o.b'
    @inbounds for i in 1:size(o.A, 1)
        for j in 1:i
            o.C[j, i] = o.A[j, i] - o.b[i] * o.b[j]
        end

    end
    LinearAlgebra.copytri!(o.C, 'U')
    o.C
end

function Statistics.cov(o::WeightedCovMatrix; corrected = false, weight_type = :analytic)
    if corrected
        if weight_type == :analytic
            LinearAlgebra.rmul!(
                value(o), 1 / (1 - (o.W2 * nobs(o)) / (weightsum(o) ^ 2))
            )
        elseif weight_type == :frequency
            LinearAlgebra.rmul!(
                value(o), 1 / (weightsum(o) - 1) * weightsum(o)
            )
        elseif weight_type == :probability
            error("If you need this, please make a PR or open an issue")
        else
            throw(ArgumentError("weight type $weight_type not implemented"))
        end
    else
        value(o)
    end
end

function Statistics.cor(o::WeightedCovMatrix; kw...)
    cov(o; kw...)
    v = diag(o.C)
    v .= 1 ./ sqrt.(v)
    return o.C .* v .* v'
end

Base.sum(o::WeightedCovMatrix) = o.b .* (meanweight(o) * nobs(o))
Statistics.mean(o::WeightedCovMatrix) = copy(o.b)
Statistics.var(o::WeightedCovMatrix; kw...) = diag(cov(o; kw...))
Statistics.std(o::WeightedCovMatrix; kw...) = sqrt.(var(o; kw...))

Base.eltype(o::WeightedCovMatrix{T}) where T = T
Base.copy(o::WeightedCovMatrix) =
    WeightedCovMatrix(copy(o.C), copy(o.A), copy(o.b), o.W, o.W2, o.n)

Base.size(x::WeightedCovMatrix, i) = size(x.C, i)
Base.size(x::WeightedCovMatrix) = size(x.C)

function Base.convert(::Type{WeightedCovMatrix{T}}, o::WeightedCovMatrix) where T
    WeightedCovMatrix{T}(
        convert(Matrix{T}, o.C),
        convert(Matrix{T}, o.A),
        convert(Vector{T}, o.b),
        convert(T,         o.W),
        convert(T,         o.W2),
        o.n
    )
end


text:
using MIRT
using Documenter
#using DemoCards
using Literate

# examples_templates, examples_theme = cardtheme("grid")
# examples, examples_cb = makedemos("examples", examples_templates)

# generate examples using Literate
lit = joinpath(@__DIR__, "lit")
src = joinpath(@__DIR__, "src")
notebooks = joinpath(src, "notebooks")

ENV["GKS_ENCODING"] = "utf-8"

DocMeta.setdocmeta!(MIRT, :DocTestSetup, :(using MIRT); recursive=true)

execute = true # Set to true for executing notebooks and documenter
nb = false # Set to true to generate the notebooks
for (root, _, files) in walkdir(lit), file in files
    splitext(file)[2] == ".jl" || continue
    ipath = joinpath(root, file)
    opath = splitdir(replace(ipath, lit=>src))[1]
    Literate.markdown(ipath, opath, documenter = execute)
    nb && Literate.notebook(ipath, notebooks, execute = execute)
end

# Documentation structure
ismd(f) = splitext(f)[2] == ".md"
pages(folder) =
    [joinpath(folder, f) for f in readdir(joinpath(src, folder)) if ismd(f)]

@show isci = get(ENV, "CI", nothing) == "true"

format = Documenter.HTML(;
    prettyurls = isci,
#   edit_link = "master",
#   assets = [examples_theme],
)

makedocs(;
    modules = [MIRT],
    authors = "Jeff Fessler and contributors",
    sitename = "MIRT.jl",
    format,
    pages = [
        "Home" => "index.md",
        "Table of Contents" => "toc.md",
        # "Examples" => examples,
        "Examples" => pages("examples"),
        "Function References" => "reference.md",
    ],
)

# examples_cb()

if isci
    deploydocs(;
        repo = "github.com/JeffFessler/MIRT.jl.git",
        devbranch = "master",
        versions = ["stable" => "v^", "dev" => "dev"],
    )
end


text:
#contains methods and associated helper functions for forming a MetaDictionary object
#each object contains the complete metadata of a file

mutable struct MetaDictionary
  fileName::String
  metaTuples::Vector{MetaTuple}
  dictionary::Dict
  index::Dict
end

#return a metadictionary object. The actual dictionary can be passed as a type lookup file
function MetaDictionary(fileName::String, metaTuples::Vector{MetaTuple})::MetaDictionary
  NTuples = length(metaTuples)
  dictionary::Dict = Dict(metaTuples[i].fieldName => metaTuples[i].columnType for i::Int ∈ 1:NTuples)
  index::Dict = Dict(metaTuples[i].fieldName => i for i::Int ∈ 1:NTuples)

  return MetaDictionary(fileName, metaTuples, dictionary, index)
end

#createas a MetaDictionary from the metadata file
function MetaDictionary(fileName::String; metaPath::String = META_PATH,
  subTablePath::String = SUB_TABLES_PATH)::MetaDictionary

  #first open the file
  iStream = open("$metaPath\\$fileName.csv")
  metaDF::DataFrame = CSV.read(iStream, strings=:raw, categorical=false)
  close(iStream)

  NFields::Int = size(metaDF,1)

  #this is where we will put the results
  metadata::Vector{MetaTuple} = Vector{MetaTuple}(NFields)

  #these are the data
  names::Vector{Symbol} = (Symbol).(metaDF[:name])
  types::Vector{Symbol} = (Symbol).(metaDF[:type])
  subTables::Vector{Bool} = (s::String->s=="True"?true:false).(metaDF[:subTable])
  shortDescriptions::Vector{MString} = metaDF[:shortDesc]
  longDescriptions::Vector{MString} = metaDF[:longDesc]

  for i::Int ∈ 1:NFields
    metadata[i] = MetaTuple(fileName, names[i], types[i],
      subTables[i], shortDescriptions[i], longDescriptions[i])
  end

  return MetaDictionary(fileName, metadata)
end

#removes fields from a meta dictionary
function removeField!(meta::MetaDictionary, field::Symbol)::MetaDictionary
  target::Int = meta.index[field]
  meta.metaTuples = meta.metaTuples[(i::Int->i≠target).(1:end)]

  NTuples = length(meta.metaTuples)
  delete!(meta.dictionary, field)
  meta.index = Dict(meta.metaTuples[i].fieldName => i for i::Int ∈ 1:NTuples)

  return meta
end

#removes fields from a meta dictionary
function removeFields!(meta::MetaDictionary, fields::Vector{Symbol})::MetaDictionary
  targets::Vector{Int} = ((s::Symbol)->meta.index[s]).(fields) #get the indices
  meta.metaTuples = meta.metaTuples[(i::Int->i ∉ targets).(1:end)] #delete the data

  NTuples = length(meta.metaTuples) #update the dictionaries
  #(delete!).(meta.dictionary, fields)
  meta.dictionary = Dict(
    meta.metaTuples[i].fieldName => meta.metaTuples[i].columnType for i::Int ∈ 1:NTuples)
  meta.index = Dict(
    meta.metaTuples[i].fieldName => i for i::Int ∈ 1:NTuples)

  return meta
end

function removeDictionary(metadata::Vector{MetaDictionary}, fileName::String)
  targets::Vector{Int} = Vector{Int}()
  NMeta::Int = length(metadata)

  for i::Int ∈ 1:NMeta
    if metadata[i].fileName == fileName
      push!(targets, i)
    end
  end

  if length(targets) < 1
    error("Attempted to delete $fileName but $fileName not found.")
  end

  metadata = metadata[((i::Int)->i ∉ targets).(1:NMeta)]

  return metadata
end

#returns a dictionary for converting data names to meta names
meta2Index(metadata::Vector{MetaDictionary})::Dict =
  Dict(metadata[i].fileName => i for i::Int ∈ 1:length(metadata))

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using JuMP
using MathOptInterface
using GLPK
using MathOptSetDistances
const MOI = MathOptInterface
const MOD = MathOptSetDistances
model = Model(GLPK.Optimizer)
@variable(model, x>=0)
@variable(model, y)
@constraint(model, c1, 2x + y <= 4)
@constraint(model, c2, x + 2y -5 == 0)
@objective(model, Max, 3x + 3y)

f = JuMP.MOI.get(m, JuMP.MOI.ConstraintFunction(), c2)
s = JuMP.MOI.get(m, JuMP.MOI.ConstraintSet(), c2)

delete(m,c2)

optimize!(m)
termination_status(m)
objective_value(m)

z = 0
z += f.constant
for t in f.terms
    coef = t.coefficient
    indx = t.variable_index.value
    global z += coef * value(all_variables(m)[indx])
end

@show z

# for (F, S) in list_of_constraint_types(m)
#     for con in all_constraints(m, F, S)
#         f = JuMP.MOI.get(m, JuMP.MOI.ConstraintFunction(), con)
#         s = JuMP.MOI.get(m, JuMP.MOI.ConstraintSet(), con)

#         delete(m,con)
#         unregister(m,:con)
#         @show s
#     end
# end


text:
# Types -------

type User <: Owner
    id
    email
    name
    login
    blog
    bio
    location
    gravatar_id
    avatar_url

    public_repos
    owned_private_repos
    total_private_repos

    public_gists
    private_gists

    followers
    following
    collaborators

    company
    hireable

    updated_at
    created_at

    plan
    user_type
    site_admin
    disk_usage

    function User(data::Dict)
        new(get(data, "id", nothing),
            get(data, "email", nothing),
            get(data, "name", nothing),
            get(data, "login", nothing),
            get(data, "blog", nothing),
            get(data, "bio", nothing),
            get(data, "location", nothing),
            get(data, "gravatar_id", nothing),
            get(data, "avatar_url", nothing),
            get(data, "public_repos", nothing),
            get(data, "owned_private_repos", nothing),
            get(data, "total_private_repos", nothing),
            get(data, "public_gists", nothing),
            get(data, "private_gists", nothing),
            get(data, "followers", nothing),
            get(data, "following", nothing),
            get(data, "collaborators", nothing),
            get(data, "company", nothing),
            get(data, "hireable", nothing),
            get(data, "updated_at", nothing),
            get(data, "created_at", nothing),
            get(data, "plan", nothing),
            get(data, "type", nothing),
            get(data, "site_admin", nothing),
            get(data, "disk_usage", nothing))
    end
end

function show(io::IO, user::User)
    print(io, "$User - $(user.login)")

    if user.name != nothing && !isempty(user.name) && user.email != nothing && !isempty(user.email)
        print(io, " ($(user.name), $(user.email))")
    elseif user.name != nothing
        print(io, " ($(user.name))")
    elseif user.email != nothing
        print(io, " ($(user.email))")
    end

    if user.bio != nothing && !isempty(user.bio)
        print(io, "\n\"$(user.bio)\"")
    end
end

# Interface -------

function user(username; auth = AnonymousAuth(), options...)
    user(auth, username; options...)
end

function user(auth::Authorization, username; headers = Dict(), options...)
    authenticate_headers(headers, auth)
    r = get(URI(API_ENDPOINT; path = "/users/$username");
            headers = headers,
            options...)

    handle_error(r)

    User(JSON.parse(r.data))
end

function followers(user::String; auth = AnonymousAuth(), options...)
    followers(auth, user; options...)
end

function followers(user::User; auth = AnonymousAuth(), options...)
    followers(auth, user.login; options...)
end

function followers(auth::Authorization, user; headers = Dict(), options...)
    authenticate_headers(headers, auth)
    r = get(URI(API_ENDPOINT; path = "/users/$user/followers");
            headers = headers,
            options...)

    handle_error(r)

    map!(f -> User(f), JSON.parse(r.data))
end

function following(user::String; auth = AnonymousAuth(), options...)
    following(auth, user; options...)
end

function following(user::User; auth = AnonymousAuth(), options...)
    following(auth, user.login; options...)
end

function following(auth::Authorization, user; headers = Dict(), options...)
    authenticate_headers(headers, auth)
    r = get(URI(API_ENDPOINT; path = "/users/$user/following");
            headers = headers,
            options...)

    handle_error(r)

    map!(f -> User(f), JSON.parse(r.data))
end


text:
mutable struct CTPF <: TopicModel
	K::Int
	M::Int
	V::Int
	U::Int
	N::Vector{Int}
	C::Vector{Int}
	R::Vector{Int}
	corp::Corpus
	topics::VectorList{Int}
	scores::Matrix{Float64}
	libs::VectorList{Int}
	drecs::VectorList{Int}
	urecs::VectorList{Int}
	a::Float64
	b::Float64
	c::Float64
	d::Float64
	e::Float64
	f::Float64
	g::Float64
	h::Float64
	alef::Matrix{Float64}
	bet::Vector{Float64}
	gimel::VectorList{Float64}
	dalet::Vector{Float64}
	he::Matrix{Float64}
	vav::Vector{Float64}
	zayin::VectorList{Float64}
	het::Vector{Float64}
	phi::Matrix{Float64}
	xi::Matrix{Float64}
	newalef::Matrix{Float64}
	newhe::Matrix{Float64}
	elbo::Float64
	newelbo::Float64

	function CTPF(corp::Corpus, K::Integer, basemodel::Union{Void, BaseTopicModel}=nothing)
		@assert ispositive(K)
		@assert !isempty(corp)
		checkcorp(corp)

		M, V, U = size(corp)
		N = [length(doc) for doc in corp]
		C = [size(doc) for doc in corp]
		R = [length(doc.readers) for doc in corp]

		@assert ispositive(U)
		@assert isequal(collect(1:U), sort(collect(keys(corp.users))))
		libs = [Int[] for _ in 1:U]
		for u in 1:U, d in 1:M
			u in corp[d].readers && push!(libs[u], d)
		end

		a, b, c, d, e, f, g, h = fill(0.1, 8)

		if isa(basemodel, Union{AbstractLDA, AbstractCTM})
			@assert isequal(size(basemodel.beta), (K, V))
			alef = exp.(basemodel.beta)
			topics = basemodel.topics
		elseif isa(basemodel, Union{AbstractfLDA, AbstractfCTM})
			@assert isequal(size(basemodel.fbeta), (K, V))
			alef = exp.(basemodel.fbeta)
			topics = basemodel.topics
		else
			alef = exp.(rand(Dirichlet(V, 1.0), K)' - 0.5)
			topics = [collect(1:V) for _ in 1:K]
		end

		bet = ones(K)
		gimel = [ones(K) for _ in 1:M]
		dalet = ones(K)
		he = ones(K, U)
		vav = ones(K)
		zayin = [ones(K) for _ in 1:M]
		het = ones(K)
		phi = ones(K, N[1]) / K
		xi = ones(2K, R[1]) / 2K

		model = new(K, M, V, U, N, C, R, copy(corp), topics, zeros(M, U), libs, Vector[], Vector[], a, b, c, d, e, f, g, h, alef, bet, gimel, dalet, he, vav, zayin, het, phi, xi)
		fixmodel!(model, check=false)

		for d in 1:M
			model.phi = ones(K, N[d]) / K
			model.xi = ones(2K, R[d]) / 2K
			updateNewELBO!(model, d)
		end
		model.phi = ones(K, N[1]) / K
		model.xi = ones(2K, R[1]) / 2K
		updateELBO!(model)
		return model
	end
end

function Elogpya(model::CTPF, d::Int)
	x = 0
	readers, ratings = model.corp[d].readers, model.corp[d].ratings
	for (u, (re, ra)) in enumerate(zip(readers, ratings)), i in 1:model.K
		binom = Binomial(ra, model.xi[i,u])
		x += (ra * model.xi[i,u] * (digamma(model.gimel[d][i]) - log(model.dalet[i]) + digamma(model.he[i,re]) - log(model.vav[i]))
				- (model.gimel[d][i] / model.dalet[i]) * (model.he[i,re] / model.vav[i]) - sum([pdf(binom, y) * lgamma(y + 1) for y in 0:ra]))
	end
	return x
end

function Elogpyb(model::CTPF, d::Int)
	x = 0
	readers, ratings = model.corp[d].readers, model.corp[d].ratings
	for (u, (re, ra)) in enumerate(zip(readers, ratings)), i in 1:model.K
		binom = Binomial(ra, model.xi[model.K+i,u])
		x += (ra * model.xi[model.K+i,u] * (digamma(model.zayin[d][i]) - log(model.het[i]) + digamma(model.he[i,re]) - log(model.vav[i]))
				- (model.zayin[d][i] / model.het[i]) * (model.he[i,re] / model.vav[i]) - sum([pdf(binom, y) * lgamma(y + 1) for y in 0:ra]))
	end
	return x
end

function Elogpz(model::CTPF, d::Int)
	x = 0
	terms, counts = model.corp[d].terms, model.corp[d].counts
	for (n, (j, c)) in enumerate(zip(terms, counts)), i in 1:model.K
		binom = Binomial(c, model.phi[i,n])
		x += (c * model.phi[i,n] * (digamma(model.gimel[d][i]) - log(model.dalet[i]) + digamma(model.alef[i,j]) - log(model.bet[i]))
				- (model.gimel[d][i] / model.dalet[i]) * (model.alef[i,j] / model.bet[i]) - sum([pdf(binom, z) * lgamma(z + 1) for z in 0:c]))
	end
	return x
end

function Elogpbeta(model::CTPF)
	x = model.V * model.K * (model.a * log(model.b) - lgamma(model.a))
	for j in 1:model.V, i in 1:model.K
		x += (model.a - 1) * (digamma(model.alef[i,j]) - log(model.bet[i])) - model.b * model.alef[i,j] / model.bet[i]
	end
	return x
end

function Elogptheta(model::CTPF, d::Int)
	x = model.K * (model.c * log(model.d) - lgamma(model.c))
	for i in 1:model.K
		x += (model.c - 1) * (digamma(model.gimel[d][i]) - log(model.dalet[i])) - model.d * model.gimel[d][i] / model.dalet[i]
	end
	return x
end

function Elogpeta(model::CTPF)
	x = model.U * model.K * (model.e * log(model.f) - lgamma(model.e))
	for u in 1:model.U, i in 1:model.K
		x += (model.e - 1) * (digamma(model.he[i,u]) - log(model.vav[i])) - model.f * model.he[i,u] / model.vav[i]
	end
	return x
end

function Elogpepsilon(model::CTPF, d::Int)
	x = model.K * (model.g * log(model.h) - lgamma(model.g))
	for i in 1:model.K
		x += (model.g - 1) * (digamma(model.zayin[d][i]) - log(model.het[i])) - model.h * model.zayin[d][i] / model.het[i]
	end
	return x
end

function Elogqy(model::CTPF, d::Int)
	x = 0
	for (u, ra) in enumerate(model.corp[d].ratings)
		x -= entropy(Multinomial(ra, model.xi[:,u]))
	end
	return x
end

function Elogqz(model::CTPF, d::Int)
	x = 0
	for (n, c) in enumerate(model.corp[d].counts)
		x -= entropy(Multinomial(c, model.phi[:,n]))
	end
	return x
end

function Elogqbeta(model::CTPF)
	x = 0
	for j in 1:model.V, i in 1:model.K
		x -= entropy(Gamma(model.alef[i,j], 1 / model.bet[i]))
	end
	return x
end

function Elogqtheta(model::CTPF, d::Int)
	x = 0
	for i in 1:model.K
		x -= entropy(Gamma(model.gimel[d][i], 1 / model.dalet[i]))
	end
	return x
end

function Elogqeta(model::CTPF)
	x = 0
	for u in 1:model.U, i in 1:model.K
		x -= entropy(Gamma(model.he[i,u], 1 / model.vav[i]))
	end
	return x
end

function Elogqepsilon(model::CTPF, d::Int)
	x = 0
	for i in 1:model.K
		x -= entropy(Gamma(model.zayin[d][i], 1 / model.het[i]))
	end
	return x
end

function updateELBO!(model::CTPF)
	model.elbo = model.newelbo + Elogpbeta(model) + Elogpeta(model) - Elogqbeta(model) - Elogqeta(model)
	model.newelbo = 0
	return model.elbo
end

function updateNewELBO!(model::CTPF, d::Int)
	model.newelbo += (Elogpya(model, d)
					+ Elogpyb(model, d)
					+ Elogpz(model, d)
					+ Elogptheta(model, d)
					+ Elogpepsilon(model, d)
					- Elogqy(model, d)
					- Elogqz(model, d)
					- Elogqtheta(model, d)
					- Elogqepsilon(model, d))
end

function updateAlef!(model::CTPF)
	model.alef = copy(model.newalef)
	model.newalef = fill(model.a, model.K, model.V)
end

function updateNewAlef!(model::CTPF, d::Int)
	terms, counts = model.corp[d].terms, model.corp[d].counts
	model.newalef[:,terms] += model.phi .* counts'
end

function updateBet!(model::CTPF)
	model.bet = model.b + sum(model.gimel) ./ model.dalet
end

function updateGimel!(model::CTPF, d::Int)
	counts, ratings = model.corp[d].counts, model.corp[d].ratings
	model.gimel[d] = model.c + model.phi * counts + model.xi[1:model.K,:] * ratings
end

function updateDalet!(model::CTPF)
	model.dalet = model.d + vec(sum(model.alef, 2)) ./ model.bet + vec(sum(model.he, 2)) ./ model.vav
end

function updateHe!(model::CTPF)
	model.he = copy(model.newhe)
	model.newhe = fill(model.e, model.K, model.U)
end

function updateNewHe!(model::CTPF, d::Int)
	readers, ratings = model.corp[d].readers, model.corp[d].ratings
	model.newhe[:,readers] += (model.xi[1:model.K,:] + model.xi[model.K+1:end,:]) .* ratings'
end

function updateVav!(model::CTPF)
	model.vav = model.f + sum(model.gimel) ./ model.dalet + sum(model.zayin) ./ model.het
end

function updateZayin!(model::CTPF, d::Int)
	ratings = model.corp[d].ratings
	model.zayin[d] = model.g + model.xi[model.K+1:end,:] * ratings
end

function updateHet!(model::CTPF)
	model.het = model.h + vec(sum(model.he, 2)) ./ model.vav
end

function updatePhi!(model::CTPF, d::Int)
	terms = model.corp[d].terms
	model.phi = exp.(digamma.(model.gimel[d]) - log.(model.dalet) - log.(model.bet) .+ digamma.(model.alef[:,terms]))
	model.phi ./= sum(model.phi, 1)
end

function updateXi!(model::CTPF, d::Int)
	readers = model.corp[d].readers
	model.xi = vcat(exp.(digamma.(model.gimel[d]) - log.(model.dalet) - log.(model.vav) .+ digamma.(model.he[:,readers])),
					exp.(digamma.(model.zayin[d]) - log.(model.het) - log.(model.vav) .+ digamma.(model.he[:,readers])))
	model.xi ./= sum(model.xi, 1)
end

function train!(model::CTPF; iter::Int=150, tol::Real=1.0, viter::Int=10, vtol::Real=1/model.K^2, chkelbo::Int=1)
	@assert all(.!isnegative.([tol, vtol]))
	@assert all(ispositive.([iter, viter, chkelbo]))

	for k in 1:iter
		chk = (k % chkelbo == 0)
		for d in 1:model.M
			for _ in 1:viter
				oldgimel = model.gimel[d]
				updatePhi!(model, d)
				updateXi!(model, d)
				updateGimel!(model, d)
				updateZayin!(model, d)
				if norm(oldgimel - model.gimel[d]) < vtol
					break
				end
			end
			updateNewAlef!(model, d)
			updateNewHe!(model, d)
			chk && updateNewELBO!(model, d)
		end
		if checkELBO!(model, k, chk, tol)
			updateDalet!(model)
			updateHet!(model)
			updateAlef!(model)
			updateBet!(model)
			updateHe!(model)
			updateVav!(model)
			break
		end
		updateDalet!(model)
		updateHet!(model)
		updateAlef!(model)
		updateBet!(model)
		updateHe!(model)
		updateVav!(model)
	end
	updatePhi!(model, 1)
	updateXi!(model, 1)

	Ebeta = model.alef ./ model.bet
	model.topics = [reverse(sortperm(vec(Ebeta[i,:]))) for i in 1:model.K]

	Eeta = model.he ./ model.vav
	for d in 1:model.M
		Etheta = model.gimel[d] ./ model.dalet
		Eepsilon = model.zayin[d] ./ model.het
		model.scores[d,:] = sum(Eeta .* (Etheta + Eepsilon), 1)
	end

	model.drecs = Vector{Int}[]
	for d in 1:model.M
		nr = setdiff(keys(model.corp.users), model.corp[d].readers)
		push!(model.drecs, nr[reverse(sortperm(vec(model.scores[d,nr])))])
	end

	model.urecs = Vector{Int}[]
	for u in 1:model.U
		ur = filter(d -> !(u in model.corp[d].readers), collect(1:model.M))
		push!(model.urecs, ur[reverse(sortperm(model.scores[ur,u]))])
	end
	nothing
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    version

function version()
    project = joinpath(dirname(@__FILE__), "..", "..", "Project.toml")
    for line in eachline(project)
        if startswith(line, "version")
            return strip(strip(split(line, " = ")[2]), ['"'])
        end
    end
end


text:
immutable IPixel
x::Int64
y::Int64
color::Int64
end

max_stars_repo_path: Module 2/Chapter03/immutable.jl

max_issues_repo_path: Module 2/Chapter03/immutable.jl

max_forks_repo_path: Module 2/Chapter03/immutable.jl