text:
using OrdinaryDiffEq
using Trixi

# Define new structs inside a module to allow re-evaluating the file.
# This module name needs to be unique among all examples, otherwise Julia will throw warnings
# if multiple test cases using the same module name are run in the same session.
module TrixiExtensionEulerRotated

using Trixi

# initial_condition_convergence_test transformed to the rotated rectangle
struct InitialConditionSourceTermsRotated
  sin_alpha::Float64
  cos_alpha::Float64
end

function InitialConditionSourceTermsRotated(alpha)
  sin_alpha, cos_alpha = sincos(alpha)

  InitialConditionSourceTermsRotated(sin_alpha, cos_alpha)
end

function (initial_condition::InitialConditionSourceTermsRotated)(x, t, equations::CompressibleEulerEquations2D)
  sin_ = initial_condition.sin_alpha
  cos_ = initial_condition.cos_alpha

  # Rotate back to unit square and translate from [-1, 1]^2 to [0, 2]^2

  # Clockwise rotation by α and translation by 1
  # Multiply with [  cos(α)  sin(α);
  #                 -sin(α)  cos(α)]
  x1 =  cos_ * x[1] + sin_ * x[2] + 1
  x2 = -sin_ * x[1] + cos_ * x[2] + 1

  rho, rho_v1, rho_v2, rho_e = initial_condition_convergence_test(SVector(x1, x2), t, equations)

  # Rotate velocity vector counterclockwise
  # Multiply with [ cos(α)  -sin(α);
  #                 sin(α)   cos(α)]
  rho_v1_rot = cos_ * rho_v1 - sin_ * rho_v2
  rho_v2_rot = sin_ * rho_v1 + cos_ * rho_v2

  return SVector(rho, rho_v1_rot, rho_v2_rot, rho_e)
end

@inline function (source_terms::InitialConditionSourceTermsRotated)(u, x, t, equations::CompressibleEulerEquations2D)
  sin_ = source_terms.sin_alpha
  cos_ = source_terms.cos_alpha

  # Rotate back to unit square and translate from [-1, 1]^2 to [0, 2]^2

  # Clockwise rotation by α and translation by 1
  # Multiply with [  cos(α)  sin(α);
  #                 -sin(α)  cos(α)]
  x1 =  cos_ * x[1] + sin_ * x[2] + 1
  x2 = -sin_ * x[1] + cos_ * x[2] + 1

  du1, du2, du3, du4 = source_terms_convergence_test(u, SVector(x1, x2), t, equations)

  # Rotate velocity vector counterclockwise
  # Multiply with [ cos(α)  -sin(α);
  #                 sin(α)   cos(α)]
  du2_rotated = cos_ * du2 - sin_ * du3
  du3_rotated = sin_ * du2 + cos_ * du3

  return SVector(du1, du2_rotated, du3_rotated, du4)
end

end # module TrixiExtensionEulerRotated

import .TrixiExtensionEulerRotated

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations2D(1.4)

alpha = 0.1
initial_condition_source_terms = TrixiExtensionEulerRotated.InitialConditionSourceTermsRotated(alpha)
sin_ = initial_condition_source_terms.sin_alpha
cos_ = initial_condition_source_terms.cos_alpha
T = [cos_ -sin_; sin_ cos_]

solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

mapping(xi, eta) = T * SVector(xi, eta)

cells_per_dimension = (16, 16)

mesh = StructuredMesh(cells_per_dimension, mapping)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_source_terms, solver,
                                    source_terms=initial_condition_source_terms)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 2.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=1.0)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
# query_string = "INSERT INTO appuser (lastname,login)
#                 VALUES (\$1, \$2)"
function create_entity!(new_object::IEntity,
                       table_name::String,
                       columns_selection_and_mapping::Dict{Symbol,<:Any},
                       dbconn::LibPQ.Connection;
                       creator::Union{IEntity,Missing} = missing)

    # tick()

    data_type = typeof(new_object)
    orm_module = PostgresORM.get_orm(new_object)

    query_string = "INSERT INTO " * table_name

    #
    # Enrich the created object with the creator and creation time if needed
    #
    creator_property = util_get_creator_property(orm_module)
    if !ismissing(creator_property)
         # Only set the creator property if it is empty because the rest of the
         #  application may want to set it (eg. if importing data)
         if ismissing(getproperty(new_object,creator_property))
             setproperty!(new_object,
                          creator_property,
                          creator)
         # If the object has the creator property set, use it
         else
             creator = getproperty(new_object,creator_property)
         end
    end

    creation_time_property = util_get_creation_time_property(orm_module)
    if !ismissing(creation_time_property)
        # Only set the creation time if it is empty because the rest of the
        #  application may want to set it (eg. if importing data)
        if ismissing(getproperty(new_object, creation_time_property))
            setproperty!(new_object,
                         creation_time_property,
                         now(Dates.UTC))
        end
    end

    # Record the state of the attributes in the modification table
    if util_get_track_changes(orm_module)

         #
         # Create the modification entries
         #
         action_id = uuid4()

         # Get the editor id if any
         editor_id = if ismissing(creator)
           missing
         else
             util_prettify_id_values_as_str(creator)
         end

         modifs = Modification[]

         previous_state_props = Dict{Symbol,Any}()

         props_modified =
            diff_dict(previous_state_props,
                      util_get_entity_props_for_comparison(new_object, dbconn))

         for (k,v) in props_modified

             modif = Modification(entity_type =  string(data_type),
                                  entity_id = missing, # not available at this moment
                                  attrname = string(k),
                                  oldvalue = ismissing(v.old) ? missing : string(v.old),
                                  newvalue = ismissing(v.new) ? missing : string(v.new),
                                  appuser_id = editor_id,
                                  creation_time = now(Dates.UTC),
                                  action_type = CRUDType.create,
                                  action_id = action_id)
             # Store the modification in a vector for the moment, we need the ID
             #  before we can persist the modifications
             push!(modifs, modif)
             # ModificationORM.create_modification(modif)

         end

    end # ENDOF if track_changes

    props = util_get_entity_props_for_db_actions(new_object,
                                                 dbconn,
                                                 false # no need to include missing values at insertion
                                                  )

    # Prepare the vectors used for creating the query
    properties_names = collect(keys(props))
    column_names = util_getcolumns(properties_names, columns_selection_and_mapping)
    properties_values = collect(values(props))

    # Loop over the properties of the object and
    #   build the appropriate list of columns
    if length(props) > 0
        query_string *= "(" * join(column_names,",") * ")" # '(lastname,login)'
    else
        # TODO throw_error
    end

    # Add the prepared statement indexes
    query_indexes = string.(collect(1:length(properties_names)))
    query_indexes = string.("\$",query_indexes)
    query_string *= (" VALUES ("
                     * join(query_indexes,",")
                     * ") " # '(lastname,login)'
                     * "RETURNING *");

    prepared_query = LibPQ.prepare(dbconn,
                                   query_string)

# @info "query is prepared"
# laptimer()

    # @info query_string
    # @info properties_values

    query_result = execute(prepared_query,
                           properties_values
                           ;throw_error=true)

# @info "query is executed"
# laptimer()

     # Get the result as a dataframe
     data = DataFrame(query_result)

     # Close the result
     close(query_result)

 # @info "query is closed"
 # laptimer()

     result = util_clean_and_rename_query_result(data,columns_selection_and_mapping)

 # @info "columns are renamed"
 # laptimer()

     result = dataframe2vector_of_namedtuples(result)

 # @info "dataframe is transformed to vector of named tuple"
 # laptimer()

     result =
        util_convert_namedtuple_to_object.(result,data_type,
                                           true, # retrieve_complex_props
                                           dbconn)

# @info "vector of named tuple is transformed to vector of objects"
# laptimer()

     # Returns the first element only because we create only one object
     result = result[1]

     util_overwrite_props!(result, # src entity
                           new_object, # target entity
                           true # Ignore missing properties like all the
                                #    unmapped properties)
                           )

     # Enrich the modification with the entity ID and persist it to database
     if util_get_track_changes(orm_module)

         idprops = util_get_ids_props_names(result)

         # Create a string representation of the ID properties values.
         id_value_as_str = util_prettify_id_values_as_str(result)

         for m in modifs
             setproperty!(m,
                          :entity_id,
                          id_value_as_str)

             # If the attribute is the ID of the entity we set the ID as the 'newvalue'
             # if m.attrname == string(orm_module.id_property)
             #     setproperty!(m,
             #                  :newvalue,
             #                  id_value_as_str)
             # end

             create_entity!(m, dbconn)
         end
     end

# tock()

     return new_object
end

function create_entity!(new_object::IEntity,
                       orm_module::Module,
                       dbconn::LibPQ.Connection;
                       creator::Union{IEntity,Missing} = missing)
   create_entity!(new_object,
                 util_get_table_with_schema_name(orm_module),
                 util_get_columns_selection_and_mapping(orm_module),
                 dbconn;
                 creator = creator)
end

function create_entity!(new_object::IEntity,
                       dbconn::LibPQ.Connection;
                       creator::Union{IEntity,Missing} = missing)

   create_entity!(new_object,
                 PostgresORM.get_orm(new_object),
                 dbconn;
                 creator = creator)
end

"""
Persist to database using Postgresql 'COPY'
"""
function create_in_bulk_using_copy(entities::Vector{T},
                                   dbconn::LibPQ.Connection)  where T <: PostgresORM.IEntity

    if isempty(entities)
        @info "Empty vector was passed to create_in_bulk_using_copy. Do nothing"
        return(0)
    end

   # dummy = T() # Does not work because the vector may be a vector of abstract types
   dummy = typeof(entities[1])()
   dataType = typeof(dummy)
   orm = PostgresORM.get_orm(dummy)
   tableName = util_get_table_name(orm)
   schemaName = util_get_schema_name(orm)
   columns_selection_and_mapping = util_get_columns_selection_and_mapping(orm)

   # Loop through the entities and check if there is a value for the fields that
   #   are reputed missing so far. Remove the field from the missing fields if
   #   there is a value.
   missing_fields = [fieldnames(dataType)...] # Initialize with all the fields
   non_missing_fields = Symbol[]
   for missing_field in missing_fields
       for e in entities
           if !ismissing(getproperty(e, missing_field))
               push!(non_missing_fields,missing_field)
               break
           end
       end
   end
   filter!(x -> !(x in non_missing_fields), missing_fields)

   # Retrieve the names of the columns from a dummy instance
   column_names = begin

       props_dummy = util_get_entity_props_for_db_actions(
                   dummy,
                   dbconn,
                   true # Include missing values
                  ;exclude_props = missing_fields
                 )

       properties_names = collect(keys(props_dummy))
       util_getcolumns(properties_names, columns_selection_and_mapping)

   end

   result::Int64 = 0
   try

        rowStrings = imap(entities) do entity

            rowStringArr = []

            props = util_get_entity_props_for_db_actions(
                        entity,
                        dbconn,
                        true # Include missing values coz all lines in the COPY
                             #   must have teh ame number of elements
                       ;exclude_props = missing_fields
                      )

            properties_values = []
            for property_name in properties_names
                push!(properties_values,props[property_name])
            end

            # properties_values = collect(values(props))

            for value in properties_values
                if ismissing(value)
                    push!(rowStringArr,"")
                else
                    if isa(value,String)
                        value = replace(value, r"\n" => " ") # Replace line breaks by space
                        value = replace(value, r"\"" => "\"\"") # Escape the double quotes
                        value = "\"$value\"" # Put the strings in between quotes so that
                                             #   we can have the delimeter appear in the strings
                    end
                    push!(rowStringArr,value)
                end
            end

            oneRowStr = "$(join(rowStringArr,';'))\n"
            oneRowStr
            # "$(row.rechor_id),$(row.rechor_sm),$(row.journeeExploitation)\n"

        end #ENDOF do row

        # return(rowStrings)
        copyin =
             LibPQ.CopyIn("COPY $schemaName.$tableName($(join(column_names,',')))
                             FROM STDIN (FORMAT CSV, DELIMITER ';', QUOTE '\"');", rowStrings)

        resultOfExecution = execute(dbconn, copyin; throw_error=false)

        # Return the number of rows inserted
        result = LibPQ.num_affected_rows(resultOfExecution)

   catch e
      rethrow(e)
   end

   return result

end

ext: jl


text:
using Base.Test

include("MNL.jl")

### Generate data

# make it deterministic
srand(10)

m, n, k = 100, 100, 3
Theta = randn(m, k) * randn(k, n)
Theta = Theta .- (Theta * ones(n,1))

nobs = int(10*m*n)
K = 2 #int(.05*n) # size of sets
sets = Array(Array{Int,1}, nobs) # set presented
rows = Array(Int, nobs)      # customer who arrived
cols = Array(Int, nobs)      # product they picked
for iobs = 1:nobs
	it = sample(1:m)
	St = sample(1:n, K, replace=false, ordered=true)
	wv = WeightVec(Float64[exp(-Theta[it, j]) for j in St])
	jt = sample(wv)
	rows[iobs], cols[iobs], sets[iobs] = it, jt, St
end
data = MNLdata(nobs, m, n, sets, rows, cols)

### Test invariances

ThetaTest = randn(m, n)
u,s,v = svd(ThetaTest)
ThetaTest = u[:,1:k]*diagm(s[1:k])*v[:,1:k]'
UTest = factor(ThetaTest, k=k)

@test_approx_eq negloglik(ThetaTest, data) negloglik(ThetaTest + rand(m,1)*ones(1,n), data)
@test_approx_eq grad_negloglik(ThetaTest + rand(m,1)*ones(1,n), data) grad_negloglik(ThetaTest + rand(m,1)*ones(1,n), data)

@test_approx_eq negloglik(ThetaTest, data) negloglik(UTest, data)

@test_approx_eq nucnorm(ThetaTest) nucnorm(UTest)

# gradients satisfy some relationship too...
# @test_approx_eq grad_negloglik(ThetaTest, data) grad_negloglik(UTest, data)
# @test_approx_eq grad_nucnorm(ThetaTest) grad_nucnorm(UTest)

max_stars_repo_name: madeleineudell/FirstOrderOptimization


text:
# Module written by Ed Scheinerman, ers@jhu.edu
# distributed under terms of the MIT license

module ShowSet

import Base.string, Base.show

AbstractSet = Union{Set,IntSet}

function string(A::AbstractSet)
    elements = collect(A)
    try
        sort!(elements)
    end
    n = length(elements)
    output = "{"
    for k=1:n
        output = output * string(elements[k])
        if k<n
            output = output * ","
        end
    end
    output = output * "}"
    return output
end

show(io::IO, A::Set)    = print(io,string(A))
show(io::IO, A::IntSet) = print(io,string(A))

end # module ShowsSet


text:
using CSV
using DataFrames
using Plots
using JLD2
using DiscretisedFluidQueues
import Distributions
include("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/plot_utils.jl")
include("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/error_metrics.jl")

@halfwidth_plot_defaults()

plotlyjs()
function p()
ic = "point_mass"
ks_data = CSV.read((@__DIR__)*"/../data/kolmogorov_smirnov.csv",DataFrame)
ks_data.DG_lin_lim = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data/errors/ks_error.csv",DataFrame).DG_limit

ks_data_lwr = CSV.read((@__DIR__)*"/../data/kolmogorov_smirnov_ci_lwr.csv",DataFrame)
ks_data_lwr.DG_lin_lim = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data/errors/ks_error_lwr.csv",DataFrame).DG_limit

ks_data_upr = CSV.read((@__DIR__)*"/../data/kolmogorov_smirnovci_upr.csv",DataFrame)
ks_data_upr.DG_lin_lim = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data/errors/ks_error_upr.csv",DataFrame).DG_limit

lwr = [log10(ks_data[i,j]) .- log10(ks_data_lwr[i,j]) for i in 1:size(ks_data,1), j in 1:size(ks_data,2)]
upr = [log10(ks_data_upr[i,j]) .- log10(ks_data[i,j]) for i in 1:size(ks_data,1), j in 1:size(ks_data,2)]
ks_data = log10.(ks_data)
plot()
linestyles_vec = [:solid,:dash,:dashdot, :dot,:dash]
markerstyles_vec = [:false,:false,:false,:false,:false]
_names = ["DG"; "DG lim"; "Unif"; "QBD-RAP"; "DG lin lim"]
colours = [1;3;2;4;5]
for (c,col) in enumerate(names(ks_data))
    os = collect(1:2:21)
    (c==5)&&(os.+=1)
    plot!(log10.(os),ks_data[:,col], ribbon=(lwr[:,c],upr[:,c]),
        label=_names[c],linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2,color=colours[c])
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="KS error - CDF")
# plot!(legend=(0.8,0.97))
plot!()
@add_lines!(ks_data,("Unif","QBDRAP"),@__DIR__)

savefig((@__DIR__)*"/ks_error_formatted.pdf")
end
p()

function p()
ic = "point_mass"
L1_cell_probs_data = CSV.read((@__DIR__)*"/../data/L1_cell_probs.csv",DataFrame)
L1_cell_probs_data.DG_lin_lim = log10.(CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data/errors/l1_cell_probs_error.csv",DataFrame).DG_limit)

L1_cell_probs_data_lwr = CSV.read((@__DIR__)*"/../data/L1_cell_probs_ci_lwr.csv",DataFrame)
L1_cell_probs_data_lwr.DG_lin_lim = log10.(CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data/errors/l1_cell_probs_error_lwr.csv",DataFrame).DG_limit)

L1_cell_probs_data_upr = CSV.read((@__DIR__)*"/../data/L1_cell_probsci_upr.csv",DataFrame)
L1_cell_probs_data_upr.DG_lin_lim = log10.(CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data/errors/l1_cell_probs_error_upr.csv",DataFrame).DG_limit)

lwr = [(L1_cell_probs_data[i,j]) .- (L1_cell_probs_data_lwr[i,j]) for i in 1:size(L1_cell_probs_data,1), j in 1:size(L1_cell_probs_data,2)]
upr = [(L1_cell_probs_data_upr[i,j]) .- (L1_cell_probs_data[i,j]) for i in 1:size(L1_cell_probs_data,1), j in 1:size(L1_cell_probs_data,2)]

p = plot()
linestyles_vec = [:solid,:dash,:dashdot,:dot,:dash]
markerstyles_vec = [:false,:false,:false,:false,:false]
_names = ["DG"; "DG lim"; "Unif"; "QBD-RAP"; "DG lin lim"]
colours = [1;3;2;4;5]
for (c,col) in enumerate(names(L1_cell_probs_data))
    os = collect(1:2:21)
    (c==5)&&(os.+=1)
    plot!(log10.(os),L1_cell_probs_data[:,col], ribbon=(lwr[:,c],upr[:,c]),
        label=_names[c],linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2,color=colours[c])
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="Cell-wise error")
# plot!(legend=(0.7,0.9))
plot!()
@add_lines!(L1_cell_probs_data,("Unif","QBDRAP"),@__DIR__)

savefig((@__DIR__)*"/L1_cell_probs_error_formatted.pdf")
end
p()

# function p()
# l2_cdf_data = CSV.read((@__DIR__)*"/../data/L2.csv",DataFrame)
# l2_cdf_data_lwr = CSV.read((@__DIR__)*"/../data/L2_ci_lwr.csv",DataFrame)
# l2_cdf_data_upr = CSV.read((@__DIR__)*"/../data/L2ci_upr.csv",DataFrame)
#  lwr = [(l2_cdf_data[i,j]) .- (l2_cdf_data_lwr[i,j]) for i in 1:size(l2_cdf_data,1), j in 1:size(l2_cdf_data,2)]
#  upr = [(l2_cdf_data_upr[i,j]) .- (l2_cdf_data[i,j]) for i in 1:size(l2_cdf_data,1), j in 1:size(l2_cdf_data,2)]
# plot()
# _names = ["DG"; "DG lim"; "Unif"; "QBD-RAP"]
# markerstyles_vec = [:false,:false,:false,:false]
# colours = [1;3;2;4]
# linestyles_vec = [:solid,:dash,:dashdot, :dot]
# for (c,col) in enumerate(names(l2_cdf_data))
#     plot!(log10.(1:2:21),l2_cdf_data[:,col],ribbon=(lwr[:,c],upr[:,c]),
#         label=_names[c],linestyle=linestyles_vec[c],
#         marker=markerstyles_vec[c],
#         linewidth=2,color=colours[c])
# end
# plot!(xlabel="Dimension")
# plot!(ylabel="Error"); error_ticks!(plot!())
# plot!(title="L² error - CDF")
# plot!(legend=:outerbottomright)
# plot!()
# @add_lines!(l2_cdf_data,("Unif","QBDRAP"),@__DIR__)

# savefig((@__DIR__)*"/l2_cdf_error_formatted.pdf")
# end
# p()

function p()
ic = "point_mass"

l1_cdf_data = CSV.read((@__DIR__)*"/../data/L1.csv",DataFrame)
l1_cdf_data.DG_lin_lim = log10.(CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data/errors/l1_error.csv",DataFrame).DG_limit)

l1_cdf_data_lwr = CSV.read((@__DIR__)*"/../data/L1_ci_lwr.csv",DataFrame)
l1_cdf_data_lwr.DG_lin_lim = log10.(CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data/errors/l1_error_lwr.csv",DataFrame).DG_limit)

l1_cdf_data_upr = CSV.read((@__DIR__)*"/../data/L1ci_upr.csv",DataFrame)
l1_cdf_data_upr.DG_lin_lim = log10.(CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data/errors/l1_error_upr.csv",DataFrame).DG_limit)

lwr = [(l1_cdf_data[i,j]) .- (l1_cdf_data_lwr[i,j]) for i in 1:size(l1_cdf_data,1), j in 1:size(l1_cdf_data,2)]
upr = [(l1_cdf_data_upr[i,j]) .- (l1_cdf_data[i,j]) for i in 1:size(l1_cdf_data,1), j in 1:size(l1_cdf_data,2)]
 markerstyles_vec = [:false,:false,:false,:false,:false]

@halfwidth_plot_defaults()
plot()
_names = ["DG"; "DG lim"; "Unif"; "QBD-RAP"; "DG lin lim"]
colours = [1;3;2;4;5]
linestyles_vec = [:solid,:dash,:dashdot,:dot,:dash]
for (c,col) in enumerate(names(l1_cdf_data))
    os = collect(1:2:21)
    (c==5)&&(os.+=1)
    plot!(log10.(os),l1_cdf_data[:,col],ribbon=(lwr[:,c],upr[:,c]),
        label=_names[c],linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2,color=colours[c])
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="L¹ error - CDF")
# plot!(legend=(0.215,0.54))
# plot!(ylims=ylims(plot!()).+(-0.175,00),xlims=xlims(plot!()).+(-0.0375,0))
@add_lines!(l1_cdf_data,("Unif","QBDRAP"),@__DIR__)

savefig((@__DIR__)*"/l1_cdf_error_formatted.pdf")
end
p()
# @fullwidth_plot_defaults()
# dir = @__DIR__()
# os = 1:2:7
# plot(layout=(4,length(os)))
# linestyles_vec = [:solid,:solid,:solid,:solid]
# markerstyles_vec = [:cross,:diamond,:circle]
# xs = range(eps(),10-eps(),length=200)
# for (c,o) in enumerate(os)
#     model = BoundedFluidQueue(fill(0.0,2,2),[1.0; 0.0],10.0)

#     mesh = DGMesh(0.0:10.0,o)
#     dq = DiscretisedFluidQueue(model,mesh)
#     jldopen(dir*"/../../../coeffs_matrix.jld2") do f
#         coeffs = f["coeff_matrix"][fun_number,c,1]
#         dt = SFMDistribution(coeffs,dq)
#         plot!(xs,pdf(dt).(xs,1),subplot=c,label=(c==1 ? "DG" : false),
#             linestyle=linestyles_vec[1],ylims = (-0.1,1.1),  yticks=false,
#             # marker=markerstyles_vec[1],
#             linewidth=2,xticks=false)
#     end
#     plot!(subplot=c, yticks=(c==1 ? (0.1:0.4:1.2) : false),
#         ylabel=(c==1 ? "DG" : false), grid=false,
#         title = string("Dimension ",o))

#     mesh = DGMesh(0.0:(1.0./o):10.0,1)
#     dq = DiscretisedFluidQueue(model,mesh)
#     jldopen(dir*"/../../../coeffs_matrix.jld2") do f
#         coeffs = f["coeff_matrix"][fun_number,c,2]
#         dt = SFMDistribution(coeffs,dq)
#         plot!(xs,pdf(dt).(xs,1),subplot=c+4,label=(c==1 ? "unif." : false),
#             linestyle=linestyles_vec[2],ylims = (-0.1,1.1),  yticks=false,
#             # marker=markerstyles_vec[2],
#             linewidth=2,xticks=false)
#     end
#     (c==1)&&plot!(subplot=c+4, yticks=(c==1 ? (0.1:0.4:1.2) : false),
#         ylabel=(c==1 ? "unif." : false), grid=false)

#     mesh = DGMesh(0.0:10.0,o)
#     dq = DiscretisedFluidQueue(model,mesh)
#     jldopen(dir*"/../../../coeffs_matrix.jld2") do f
#         coeffs = f["coeff_matrix"][fun_number,c,3]
#         dt = SFMDistribution(coeffs,dq)
#         plot!(xs,pdf(dt).(xs,1),subplot=c+8,label=(c==1 ? "DG (limiter)" : false),
#             linestyle=linestyles_vec[3],ylims = (-0.1,1.1),  yticks=false,
#             # marker=markerstyles_vec[3],
#             linewidth=2,xticks=false)
#     end
#     (c==1)&&plot!(subplot=c+8, yticks=(c==1 ? (0.1:0.4:1.2) : false),
#         ylims = (-0.1,1.1),
#         ylabel=(c==1 ? "DG-limit" : false), grid=false)

#     mesh = FRAPMesh(0.0:10.0,o)
#     dq = DiscretisedFluidQueue(model,mesh)
#     jldopen(dir*"/../../../coeffs_matrix.jld2") do f
#         coeffs = f["coeff_matrix"][fun_number,c,4]
#         dt = SFMDistribution(coeffs,dq)
#         plot!(xs,pdf(dt).(xs,1),subplot=c+12,label=(c==1 ? "QBD-RAP" : false),
#             linestyle=linestyles_vec[3],ylims = (-0.1,1.1), yticks=false,
#             # marker=markerstyles_vec[3],
#             linewidth=2,xticks=0:2:10, grid=false)
#     end
#     (c==1)&&plot!(subplot=c+12, yticks=(c==1 ? (0.1:0.4:1.2) : false),
#         ylims = (-0.1,1.1),
#         ylabel=(c==1 ? "QBD-RAP" : false), grid=false)
#     # plot!(xs,x->2*x*(x<0.5) + Float64(x>=0.5),subplot=c,label=(c==1 ? "truth" : false),
#     #     linestyle=:dot,
#     #     # marker=markerstyles_vec[3],
#     #     color=:grey,
#     #     linewidth=2)

#     # plot!(subplot=c, yticks=(c==1 ? (0.1:0.3:1.2) : false),
#     #     ylims = (-0.1,1.1), xticks = 0.0:3:10,
#     #     ylabel=(c==1 ? "probability" : false), grid=false,
#     #     title = string("Dimension ",o),
#     #     legend=:outertopright)
#     plot!(legend=false)
# end
# plot!()
# savefig((@__DIR__)*"/pdf_formatted.pdf")

# @halfwidth_plot_defaults()
# L1_cell_probs_errors = DataFrame(DG = Float64[], Order_1 = Float64[], DG_limiter = Float64[], QBDRAP = Float64[])
# i=1; t=4
# F4(z) = exp(-1*t)*(exp(z)-1)*(z<=t) + (z>t)*(1-exp(-1*t))
# truth = F4.(1:10)-F4.(0:9)
# truth_pm = zeros(3)
# for (c_o,o) in enumerate(1:2:21)
#     row = zeros(4)
#     for (c_m,m) in enumerate(names(L1_cell_probs_errors))
#         jldopen(dir*"/../../../coeffs_matrix.jld2") do f
#             model = BoundedFluidQueue(fill(0.0,2,2),[1.0; 0.0],10.0)

#             mesh = (c_m<4) ? DGMesh(0.0:10.0,o) : FRAPMesh(0.0:10.0,o)
#             dq = DiscretisedFluidQueue(model,mesh)

#             coeffs = f["coeff_matrix"][fun_number,c_o,c_m]
#             dt = SFMDistribution(coeffs,dq)

#             x_vals = range(0.5,9.5;length=10)

#             approx = sum(cell_probs(dt).(x_vals,(1:2)'),dims=2)
#             approx_pm = [coeffs[1:N₋(model)];coeffs[end-N₊(model)+1:end]]

#             (o==3)&&display(approx)
#             row[c_m] = log10(sum(abs.(approx-truth)) + sum(abs.(approx_pm-truth_pm)))
#         end
#     end
#     push!(L1_cell_probs_errors,row)
# end

# file = (@__DIR__)*"/../data/meshs_l1_cell_"*"func_count_$(fun_number)"
# CSV.write(file*".csv",L1_cell_probs_errors)
# q = plot()
# linestyles_vec = [:solid,:dash,:dashdot, :dot]
# markerstyles_vec = [:cross,:diamond,:circle,:dot]
# for (c,reconstruction) in enumerate(names(L1_cell_probs_errors))
#     plot!(log10.(1:2:21),L1_cell_probs_errors[:,reconstruction],label=reconstruction,
#         linestyle=linestyles_vec[c],
#         marker=markerstyles_vec[c],
#         linewidth=2)
# end
# plot!(xlabel="Dimension", ylabel="log₁₀ Error",
#     title="Error between cell masses",
#     legend=:outerbottomright)
# savefig((@__DIR__)*"/L1_cell_probs.pdf")


text:
push!(LOAD_PATH, pwd()*"/../../src/")
using Distributions
using PyPlot
using BVB

x = linspace(-0.1, 1.1, 500)
init_approx_d = Normal(0, 4)

lower = 0.
upper = 1.
max_iter = 100

# example 1
@show true_d = Beta(2.0, 2.0)
ll(x) = logpdf(true_d, x)

# include("run_laplacian_gradient_boosting.jl")
include("run_laplacian_gradient_boosting_exclusive.jl")
show()

# example 2
@show true_d = Beta(2.0, 5.0)
ll(x) = logpdf(true_d, x)

# include("run_laplacian_gradient_boosting.jl")
include("run_laplacian_gradient_boosting_exclusive.jl")
show()

# example 3
@show true_d = Beta(0.6, 0.5)
ll(x) = logpdf(true_d, x)

# include("run_laplacian_gradient_boosting.jl")
include("run_laplacian_gradient_boosting_exclusive.jl")
show()

println("Done.")
readline(STDIN)

max_stars_repo_name: richardkwo/boosting-VB

max_issues_repo_name: richardkwo/boosting-VB

max_forks_repo_name: richardkwo/boosting-VB


text:
immutable MultivariateNormal <: ContinuousMultivariateDistribution
    mean::Vector{Float64}
    covchol::Cholesky{Float64}
    function MultivariateNormal(m::Vector{Float64}, c::Cholesky{Float64})
        if length(m) == size(c, 1) == size(c, 2)
            new(m, c)
        else
            error("Dimensions of mean vector and covariance matrix do not match")
        end
    end
end

function MultivariateNormal(mean::Vector{Float64}, cov::Matrix{Float64})
    MultivariateNormal(mean, cholfact(cov))
end
function MultivariateNormal(mean::Vector{Float64})
    MultivariateNormal(mean, eye(length(mean)))
end
function MultivariateNormal(cov::Matrix{Float64})
    MultivariateNormal(zeros(size(cov, 1)), cov)
end
MultivariateNormal() = MultivariateNormal(zeros(2), eye(2))

function cdf{T <: Real}(d::MultivariateNormal, x::Vector{T})
    k = length(d.mean)
    if k > 3
        error("Dimension larger than three is not supported yet")
    end
    stddev = sqrt(diag(var(d)))
    z = (x - d.mean) ./ stddev
    C = diagmm(d.covchol[:U], 1.0 / stddev)
    C = C'C
    if k == 3
        return tvtcdf(0, z, C[[2, 3, 6]])
    elseif k == 2
        return bvtcdf(0, z[1], z[2], C[2])
    else
        return cdf(Normal(), z[1])
    end
end

function entropy(d::MultivariateNormal)
    S = cov(d)
    n = size(S, 1)
    for j in 1:n
        for i in 1:n
            S[i, j] *= 2.0 * pi * e
        end
    end
    return 0.5 * logdet(S)
end

mean(d::MultivariateNormal) = d.mean

function mgf(d::MultivariateNormal, t::AbstractVector)
    m, S = d.mean, var(d)
    return exp(dot(t, m) + 0.5 * dot(t, S * t))
end

function cf(d::MultivariateNormal, t::AbstractVector)
    m, S = d.mean, var(d)
    return exp(im * dot(t, m) - 0.5 * dot(t, S * t))
end

pdf{T<:Real}(d::MultivariateNormal, x::Vector{T}) = exp(logpdf(d, x))

function logpdf{T <: Real}(d::MultivariateNormal, x::Vector{T})
    k = length(d.mean)
    u = x - d.mean
    chol_ldiv!(d.covchol, u)
    return -0.5 * (k * log(2.0 * pi) + logdet(d.covchol) + dot(u, u))
end

function logpdf!{T <: Real}(r::AbstractVector,
                            d::MultivariateNormal,
                            x::Matrix{T})
    mu::Vector{Float64} = d.mean
    k = length(mu)
    if size(x, 1) != k
        throw(ArgumentError("The dimension of x is inconsistent with d."))
    end
    n = size(x, 2)
    u = Array(Float64, k, n)
    for j = 1:n # u[:,j] = x[:,j] - mu
        for i = 1:k
            u[i, j] = x[i, j] - mu[i]
        end
    end
    chol_ldiv!(d.covchol, u)
    c::Float64 = -0.5 * (k * log(2.0 * pi) + logdet(d.covchol))
    for j = 1:n
        dot_uj = 0.0
        for i = 1:k
            dot_uj += u[i, j] * u[i, j]
        end
        r[j] = c - 0.5 * dot_uj
    end
end

function logpdf{T <: Real}(d::MultivariateNormal, x::Matrix{T})
    r = Array(Float64, size(x, 2))
    logpdf!(r, d, x)
    return r
end

function rand!(d::MultivariateNormal, x::Vector, tmp::Vector)
    randn!(tmp)
    At_mul_B(x, d.covchol[:U], tmp)
    for i in 1:length(x)
        x[i] += d.mean[i]
    end
    return x
end

function rand!(d::MultivariateNormal, x::Vector)
    tmp = randn(length(x))
    return rand!(d, x, tmp)
end

function rand(d::MultivariateNormal)
    x = Array(Float64, length(d.mean))
    tmp = Array(Float64, length(d.mean))
    return rand!(d, x, tmp)
end

function rand!(d::MultivariateNormal, X::Matrix)
    k = length(d.mean)
    tmp = Array(Float64, k)
    m, n = size(X)
    if m != k
        throw(ArgumentError("Wrong dimensions"))
    end
    randn!(X)
    for i in 1:n
        # This makes a copy of X[:, i]
        At_mul_B(tmp, d.covchol[:U], X[:, i])
        for dim in 1:m
            X[dim, i] = tmp[dim] + d.mean[dim]
        end
    end
    return X
end

function var(d::MultivariateNormal)
    U = d.covchol[:U]
    return U'U
end

function chol_ldiv!(chol::Cholesky{Float64}, u::VecOrMat{Float64})
    Base.LinAlg.LAPACK.trtrs!('U', 'T', 'N', chol.UL, u)
end


text:
@testset "binarymaximum.jl" begin

    @testset "BinaryMaximumBC()" begin
        trailer = SeaPearl.Trailer()

        x = SeaPearl.IntVar(1, 3, "x", trailer)
        y = SeaPearl.IntVar(1, 4, "y", trailer)
        z = SeaPearl.IntVar(1, 12, "z", trailer)

        constraint = SeaPearl.BinaryMaximumBC(x, y, z, trailer)

        @test constraint in x.onDomainChange
        @test constraint in y.onDomainChange
        @test constraint in z.onDomainChange
        @test constraint.active.value
    end

    @testset "propagate!(::BinaryMaximumBC)" begin

        @testset "x too small" begin

            trailer = SeaPearl.Trailer()

            x = SeaPearl.IntVar(1, 3, "x", trailer)
            y = SeaPearl.IntVar(4, 6, "y", trailer)
            z = SeaPearl.IntVar(5, 12, "z", trailer)

            constraint = SeaPearl.BinaryMaximumBC(x, y, z, trailer)
            toPropagate = Set{SeaPearl.Constraint}()
            prunedDomains = SeaPearl.CPModification()

            @test !SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        end

        @testset "z too big" begin

            trailer = SeaPearl.Trailer()

            x = SeaPearl.IntVar(2, 4, "x", trailer)
            y = SeaPearl.IntVar(1, 4, "y", trailer)
            z = SeaPearl.IntVar(8, 12, "z", trailer)

            constraint = SeaPearl.BinaryMaximumBC(x, y, z, trailer)
            toPropagate = Set{SeaPearl.Constraint}()
            prunedDomains = SeaPearl.CPModification()

            @test !SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        end

        @testset "x too big" begin

            trailer = SeaPearl.Trailer()

            x = SeaPearl.IntVar(15, 18, "x", trailer)
            y = SeaPearl.IntVar(3, 6, "y", trailer)
            z = SeaPearl.IntVar(4, 12, "z", trailer)

            constraint = SeaPearl.BinaryMaximumBC(x, y, z, trailer)
            toPropagate = Set{SeaPearl.Constraint}()
            prunedDomains = SeaPearl.CPModification()

            @test !SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        end

        @testset "Prune x below" begin

            trailer = SeaPearl.Trailer()

            x = SeaPearl.IntVar(1, 6, "x", trailer)
            y = SeaPearl.IntVar(3, 4, "y", trailer)
            z = SeaPearl.IntVar(4, 12, "z", trailer)

            constraint = SeaPearl.BinaryMaximumBC(x, y, z, trailer)
            toPropagate = Set{SeaPearl.Constraint}()
            prunedDomains = SeaPearl.CPModification()

            @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)
            @test prunedDomains == SeaPearl.CPModification("x" => [1, 2, 3])
            @test constraint.active.value

        end

        @testset "Prune x above" begin

            trailer = SeaPearl.Trailer()

            x = SeaPearl.IntVar(5, 15, "x", trailer)
            y = SeaPearl.IntVar(3, 8, "y", trailer)
            z = SeaPearl.IntVar(4, 12, "z", trailer)

            constraint = SeaPearl.BinaryMaximumBC(x, y, z, trailer)
            toPropagate = Set{SeaPearl.Constraint}()
            prunedDomains = SeaPearl.CPModification()

            @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)
            @test prunedDomains == SeaPearl.CPModification("x" => [13, 14, 15])
            @test constraint.active.value

        end

        @testset "Prune y" begin

            trailer = SeaPearl.Trailer()

            x = SeaPearl.IntVar(4, 5, "x", trailer)
            y = SeaPearl.IntVar(4, 7, "y", trailer)
            z = SeaPearl.IntVar(1, 3, "z", trailer)

            constraint = SeaPearl.BinaryMaximumBC(x, y, z, trailer)
            toPropagate = Set{SeaPearl.Constraint}()
            prunedDomains = SeaPearl.CPModification()

            @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)
            @test prunedDomains == SeaPearl.CPModification("y" => [6, 7])
            @test constraint.active.value

        end

        @testset "Prune z" begin

            trailer = SeaPearl.Trailer()

            x = SeaPearl.IntVar(7, 10, "x", trailer)
            y = SeaPearl.IntVar(1, 6, "y", trailer)
            z = SeaPearl.IntVar(5, 12, "z", trailer)

            constraint = SeaPearl.BinaryMaximumBC(x, y, z, trailer)
            toPropagate = Set{SeaPearl.Constraint}()
            prunedDomains = SeaPearl.CPModification()

            @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)
            @test prunedDomains == SeaPearl.CPModification("z" => [5, 6, 11, 12])
            @test constraint.active.value

        end

        @testset "Prune x & y" begin

            trailer = SeaPearl.Trailer()

            x = SeaPearl.IntVar(7, 10, "x", trailer)
            y = SeaPearl.IntVar(4, 8, "y", trailer)
            z = SeaPearl.IntVar(1, 4, "z", trailer)

            constraint = SeaPearl.BinaryMaximumBC(x, y, z, trailer)
            toPropagate = Set{SeaPearl.Constraint}()
            prunedDomains = SeaPearl.CPModification()

            @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)
            @test prunedDomains == SeaPearl.CPModification("x" => [9, 10], "y" => [4, 5, 6])
            @test constraint.active.value

        end

        @testset "Assign x and deactivate" begin

            trailer = SeaPearl.Trailer()

            x = SeaPearl.IntVar(7, 10, "x", trailer)
            y = SeaPearl.IntVar(8, 8, "y", trailer)
            z = SeaPearl.IntVar(1, 4, "z", trailer)

            constraint = SeaPearl.BinaryMaximumBC(x, y, z, trailer)
            toPropagate = Set{SeaPearl.Constraint}()
            prunedDomains = SeaPearl.CPModification()

            @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)
            @test SeaPearl.assignedValue(x) == 8
            @test prunedDomains == SeaPearl.CPModification("x" => [7, 9, 10])
            @test !constraint.active.value

        end

        @testset "Assign x and deactivate" begin

            trailer = SeaPearl.Trailer()

            x = SeaPearl.IntVar(1, 10, "x", trailer)
            y = SeaPearl.IntVar(3, 3, "y", trailer)
            z = SeaPearl.IntVar(4, 4, "z", trailer)

            constraint = SeaPearl.BinaryMaximumBC(x, y, z, trailer)
            toPropagate = Set{SeaPearl.Constraint}()
            prunedDomains = SeaPearl.CPModification()

            @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)
            @test SeaPearl.assignedValue(x) == 4
            @test prunedDomains == SeaPearl.CPModification("x" => [1, 2, 3, 5, 6, 7, 8, 9, 10])
            @test !constraint.active.value

        end

    end
end


text:
using BenchmarkTools, SymbolicUtils
using SymbolicUtils: isnumber

using Random

SUITE = BenchmarkGroup()

@syms a b c
let r = @rule(~x => ~x), rs = RuleSet([r]),
    acr = @rule(~x::isnumber + ~y => ~y)

    overhead = SUITE["overhead"]  = BenchmarkGroup()
    overhead["rule"]  = BenchmarkGroup()

    overhead["rule"]["noop:Int"]  = @benchmarkable $r(1)
    overhead["rule"]["noop:Sym"]  = @benchmarkable $r($a)
    overhead["rule"]["noop:Term"] = @benchmarkable $r($(a+2))

    overhead["acrule"]  = BenchmarkGroup()
    overhead["acrule"]["noop:Int"]  = @benchmarkable $acr(1)
    overhead["acrule"]["noop:Sym"]  = @benchmarkable $acr($a)
    overhead["acrule"]["a+2"] = @benchmarkable $acr($(a+2))
    overhead["acrule"]["a+b"] = @benchmarkable $acr($(a+b))
    overhead["acrule"]["a+2+b"] = @benchmarkable $acr($(a+2+b))

    overhead["ruleset"]  = BenchmarkGroup()
    overhead["ruleset"]["noop:Int"]  = @benchmarkable $rs(1)
    overhead["ruleset"]["noop:Sym"]  = @benchmarkable $rs($a)
    overhead["ruleset"]["noop:Term"] = @benchmarkable $rs($(a+2))

    overhead["simplify"]  = BenchmarkGroup()
    overhead["simplify"]["noop:Int"]  = @benchmarkable simplify(1)
    overhead["simplify"]["noop:Sym"]  = @benchmarkable simplify($a)
    overhead["simplify"]["noop:Term"] = @benchmarkable simplify($(a+2))

    overhead["simplify_no_fixp"]  = BenchmarkGroup()
    overhead["simplify_no_fixp"]["noop:Int"]  = @benchmarkable simplify(1, fixpoint=false)
    overhead["simplify_no_fixp"]["noop:Sym"]  = @benchmarkable simplify($a, fixpoint=false)
    overhead["simplify_no_fixp"]["noop:Term"] = @benchmarkable simplify($(a+2), fixpoint=false)
end


text:
# ------------ Linked map ------------ #

mutable struct LinkedMap{K} <: AddressMap{K}
    pointer::K
    next::AddressMap{<:K}
    LinkedMap{K}() where K = new{K}()
    LinkedMap(pointer::K) where K = new{K}(pointer)
    LinkedMap{K}(next::AddressMap{<:K}) where K = new{K}(next)
end


text:
s = """
@doc "this is a test" -> f() = 0
"""
msgs = lintstr(s)
@test isempty(msgs)

s = "@doc \"\"\"this is a test\"\"\" -> f() = 0"
msgs = lintstr(s)
@test isempty(msgs)

s = """
@doc "this is a test"
f() = 0
"""
msgs = lintstr(s)
@test msgs[1].code == :W443
@test contains(msgs[1].message, "did you forget an -> after @doc or make it inline?")

s = """
@doc "this is a test" f() = 0
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
@doc \"\"\"
this is a test
\"\"\" f() = 0
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
@doc f
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
func(v) = v

\"\"\"
Documentation
\"\"\"
func(v)
"""
msgs = lintstr(s)
@test isempty(msgs)

# TODO: should we warn this documentation?
"""
type TestType
    "doc"
    a
end
"""
msgs = lintstr(s)
@test isempty(msgs)


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

intersectionat(α::Float64) = Intersection(α, [0.7071067811865475, 0.7071067811865475, 0.0] * α, [0.0, 0.0, 0.0], 0.0, 0.0, NullInterface())

const greenwavelength = Unitful.u"550nm"
const redwavelength = Unitful.u"650nm"
const bluewavelength = Unitful.u"450nm"
const roomtemperature = 20 * Unitful.u"°C"

transmissivethingrating(period, orders) = ThinGratingInterface(SVector(0.0, 1.0, 0.0), period, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, minorder = -orders, maxorder = orders, reflectance = [0.0 for _ in (-orders):orders], transmission = [0.2 for _ in (-orders):orders])
reflectivethingrating(period, orders) = ThinGratingInterface(SVector(0.0, 1.0, 0.0), period, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, minorder = -orders, maxorder = orders, transmission = [0.0 for _ in (-orders):orders], reflectance = [0.2 for _ in (-orders):orders])

function multiHOE()
    rect = Rectangle(5.0, 5.0, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, 0.0))
    inbeam = SVector(0.0, 0.0, -1.0), CollimatedBeam

    int1 = HologramInterface(SVector(-5.0, 0.0, -20.0), ConvergingBeam, SVector(0.0, -1.0, -1.0), CollimatedBeam, 0.55, 100.0, OpticSim.GlassCat.Air, OpticSim.Examples.Examples_N_BK7, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, 0.05, false)
    int2 = HologramInterface(SVector(5.0, 0.0, -20.0), ConvergingBeam, SVector(0.0, 1.0, -1.0), CollimatedBeam, 0.55, 100.0, OpticSim.GlassCat.Air, OpticSim.Examples.Examples_N_BK7, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, 0.05, false)

    mint = MultiHologramInterface(int1, int2)
    obj = MultiHologramSurface(rect, mint)
    sys = CSGOpticalSystem(LensAssembly(obj), Rectangle(50.0, 50.0, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, -20.0), interface = opaqueinterface()))
    return sys
end


text:
function max_dt(u, t, mesh::TreeMesh{3},
                constant_speed::Val{false}, equations, dg::DG, cache)
  # to avoid a division by zero if the speed vanishes everywhere,
  # e.g. for steady-state linear advection
  max_scaled_speed = nextfloat(zero(t))

  for element in eachelement(dg, cache)
    max_λ1 = max_λ2 = max_λ3 = zero(max_scaled_speed)
    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
      u_node = get_node_vars(u, equations, dg, i, j, k, element)
      λ1, λ2, λ3 = max_abs_speeds(u_node, equations)
      max_λ1 = max(max_λ1, λ1)
      max_λ2 = max(max_λ2, λ2)
      max_λ3 = max(max_λ3, λ3)
    end
    inv_jacobian = cache.elements.inverse_jacobian[element]
    max_scaled_speed = max(max_scaled_speed, inv_jacobian * (max_λ1 + max_λ2 + max_λ3))
  end

  return 2 / (nnodes(dg) * max_scaled_speed)
end

function max_dt(u, t, mesh::TreeMesh{3},
                constant_speed::Val{true}, equations, dg::DG, cache)
  # to avoid a division by zero if the speed vanishes everywhere,
  # e.g. for steady-state linear advection
  max_scaled_speed = nextfloat(zero(t))

  for element in eachelement(dg, cache)
    max_λ1, max_λ2, max_λ3 = max_abs_speeds(equations)
    inv_jacobian = cache.elements.inverse_jacobian[element]
    max_scaled_speed = max(max_scaled_speed, inv_jacobian * (max_λ1 + max_λ2 + max_λ3))
  end

  return 2 / (nnodes(dg) * max_scaled_speed)
end

function max_dt(u, t, mesh::CurvedMesh{3},
                constant_speed::Val{false}, equations, dg::DG, cache)
  # to avoid a division by zero if the speed vanishes everywhere,
  # e.g. for steady-state linear advection
  max_scaled_speed = nextfloat(zero(t))

  @unpack contravariant_vectors = cache.elements

  for element in eachelement(dg, cache)
    max_λ1 = max_λ2 = max_λ3 = zero(max_scaled_speed)
    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
      u_node = get_node_vars(u, equations, dg, i, j, k, element)
      λ1, λ2, λ3 = max_abs_speeds(u_node, equations)

      Ja11, Ja12, Ja13 = get_contravariant_vector(1, contravariant_vectors, i, j, k, element)
      λ1_transformed   = abs(Ja11 * λ1 + Ja12 * λ2 + Ja13 * λ3)
      Ja21, Ja22, Ja23 = get_contravariant_vector(2, contravariant_vectors, i, j, k, element)
      λ2_transformed   = abs(Ja21 * λ1 + Ja22 * λ2 + Ja23 * λ3)
      Ja31, Ja32, Ja33 = get_contravariant_vector(3, contravariant_vectors, i, j, k, element)
      λ3_transformed   = abs(Ja31 * λ1 + Ja32 * λ2 + Ja33 * λ3)

      inv_jacobian = abs(cache.elements.inverse_jacobian[i, j, k, element])

      max_λ1 = max(max_λ1, inv_jacobian * λ1_transformed)
      max_λ2 = max(max_λ2, inv_jacobian * λ2_transformed)
      max_λ3 = max(max_λ3, inv_jacobian * λ3_transformed)
    end

    max_scaled_speed = max(max_scaled_speed, max_λ1 + max_λ2 + max_λ3)
  end

  return 2 / (nnodes(dg) * max_scaled_speed)
end

function max_dt(u, t, mesh::CurvedMesh{3},
                constant_speed::Val{true}, equations, dg::DG, cache)
  # to avoid a division by zero if the speed vanishes everywhere,
  # e.g. for steady-state linear advection
  max_scaled_speed = nextfloat(zero(t))

  @unpack contravariant_vectors = cache.elements

  max_λ1, max_λ2, max_λ3 = max_abs_speeds(equations)

  for element in eachelement(dg, cache)
    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
      Ja11, Ja12, Ja13 = get_contravariant_vector(1, contravariant_vectors, i, j, k, element)
      λ1_transformed   = abs(Ja11 * max_λ1 + Ja12 * max_λ2 + Ja13 * max_λ3)
      Ja21, Ja22, Ja23 = get_contravariant_vector(2, contravariant_vectors, i, j, k, element)
      λ2_transformed   = abs(Ja21 * max_λ1 + Ja22 * max_λ2 + Ja23 * max_λ3)
      Ja31, Ja32, Ja33 = get_contravariant_vector(3, contravariant_vectors, i, j, k, element)
      λ3_transformed   = abs(Ja31 * max_λ1 + Ja32 * max_λ2 + Ja33 * max_λ3)

      inv_jacobian = abs(cache.elements.inverse_jacobian[i, j, k, element])

      max_scaled_speed = max(max_scaled_speed,
                             inv_jacobian * (λ1_transformed + λ2_transformed + λ3_transformed))
    end
  end

  return 2 / (nnodes(dg) * max_scaled_speed)
end


text:
using GLMakie
# LAxis -> Axis
# LText -> Label
# LColorbar -> ColorBar
# LSlider -> Slider
# LMenu -> Menu
# AbstractPlotting. -> GLMakie.AbstractPlotting
# scene.children[1][Axis] -> scene.children[1][OldAxis]

aspect = (1, 1, 32/192)
resolution = (1920, 1080)
statistics = true
title = "Field = "
bins = 300

scene, layout = layoutscene(resolution = resolution)
    lscene = layout[2:4, 2:4] = LScene(scene)
    width = round(Int, resolution[1] / 4) # make menu 1/4 of preliminary resolution

    # Create choices and nodes
    stateindex = collect(1:length(states))
    statenode = Node(stateindex[1])

     colorchoices = [:balance, :thermal, :dense, :deep, :curl, :thermometer]
    colornode = Node(colorchoices[1])

    if statistics
        llscene = layout[4,1] = Axis(scene, xlabel = @lift(statenames[$statenode] * units[$statenode]),
                         xlabelcolor = :black, ylabel = "pdf",
                         ylabelcolor = :black, xlabelsize = 40, ylabelsize = 40,
                         xticklabelsize = 25, yticklabelsize = 25,
                         xtickcolor = :black, ytickcolor = :black,
                         xticklabelcolor  = :black, yticklabelcolor = :black)
        layout[3, 1] = Label(scene, "Statistics", width = width, textsize = 50)
    end

    # x,y,z are for determining the aspect ratio of the box
    if (typeof(aspect) <: Tuple) & (length(aspect) == 3)
        x, y, z = aspect
    else
        x, y, z = size(states[1])
    end

    # Clim sliders
    upperclim_slider = Slider(scene, range = range(0, 1, length = 101), startvalue = 0.99)
    upperclim_node = upperclim_slider.value
    lowerclim_slider = Slider(scene, range = range(0, 1, length = 101), startvalue = 0.01)
    lowerclim_node = lowerclim_slider.value

    # Lift Nodes
    state = @lift(states[$statenode])
    statename = @lift(statenames[$statenode])
    clims = @lift((quantile($state[:], $lowerclim_node) , quantile($state[:], $upperclim_node)))
    cmap_rgb = @lift(to_colormap($colornode))
    titlename = @lift(title * $statename) # use padding and appropriate centering

    # Statistics
    if statistics
        histogram_node = @lift(histogram($state, bins = bins))
        xs = @lift($histogram_node[1])
        ys = @lift($histogram_node[2])
        pdf = GLMakie.AbstractPlotting.barplot!(llscene, xs, ys, color = :red,
                        strokecolor = :red,
                        strokewidth = 1)
        @lift(GLMakie.AbstractPlotting.xlims!(llscene, extrema($state)))
        @lift(GLMakie.AbstractPlotting.ylims!(llscene, extrema($histogram_node[2])))
        vlines!(llscene, @lift($clims[1]), color = :black, linewidth = width / 100)
        vlines!(llscene, @lift($clims[2]), color = :black, linewidth = width / 100)
    end

    # Volume Plot
    volume!(lscene, 0..x, 0..y, 0..z, state,
            camera = cam3d!,
            colormap = cmap_rgb,
            colorrange = clims)
    # Camera
    cam = cameracontrols(scene.children[1])
    eyeposition = Float32[2, 2, 1.3]
    lookat = Float32[0.82, 0.82, 0.1]
    # Title
    supertitle = layout[1, 2:4] = Label(scene, titlename , textsize = 50, color = :black)

    # Menus
    statemenu = Menu(scene, options = zip(statenames, stateindex))
    on(statemenu.selection) do s
        statenode[] = s
    end

    colormenu = Menu(scene, options = zip(colorchoices, colorchoices))
    on(colormenu.selection) do s
        colornode[] = s
    end
    lowerclim_string = @lift("lower clim quantile = " *  @sprintf("%0.2f", $lowerclim_node) * ", value = " * @sprintf("%0.1e", $clims[1]))
    upperclim_string = @lift("upper clim quantile = " *  @sprintf("%0.2f", $upperclim_node) * ", value = " * @sprintf("%0.1e", $clims[2]))
    # depends on makie version, vbox for old, vgrid for new
    layout[2, 1] = vgrid!(
        Label(scene, "State", width = nothing),
        statemenu,
        Label(scene, "Color", width = nothing),
        colormenu,
        Label(scene, lowerclim_string, width = nothing),
        lowerclim_slider,
        Label(scene, upperclim_string, width = nothing),
        upperclim_slider,
    )
    layout[1,1] = Label(scene, "Menu", width = width, textsize = 50)

    # Modify Axis

    axis = scene.children[1][OldAxis]
    # axis[:names][:axisnames] = ("↓ Zonal [m] ", "Meriodonal [m]↓ ", "Depth [m]↓ ")
    axis[:names][:axisnames] = ("↓", "↓ ", "↓ ")

    axis[:names][:align] = ((:left, :center), (:right, :center), (:right, :center))
    # need to adjust size of ticks first and then size of axis names
    axis[:names][:textsize] = (50.0, 50.0, 50.0)
    axis[:ticks][:textsize] = (00.0, 00.0, 00.0)
    # axis[:ticks][:ranges_labels].val # current axis labels
    xticks = collect(range(-0, aspect[1], length = 2))
    yticks = collect(range(-0, aspect[2], length = 6))
    zticks = collect(range(-0, aspect[3], length = 2))
    ticks = (xticks, yticks, zticks)
    axis[:ticks][:ranges] = ticks
    xtickslabels = [@sprintf("%0.1f", (xtick)) for xtick in xticks]
    xtickslabels[end] = "1e6"
    ytickslabels = ["", "south","", "", "north", ""]
    ztickslabels = [@sprintf("%0.1f", (xtick)) for xtick in xticks]
    labels = (xtickslabels, ytickslabels, ztickslabels)
    axis[:ticks][:labels] = labels

    display(scene)

max_stars_repo_name: sandreza/Mesoscale

max_issues_repo_name: sandreza/Mesoscale

max_forks_repo_name: sandreza/Mesoscale


text:
# KKT QP with Equality Constraints Solver

include("../auglag/auglag-core.jl")

using SparseArrays

@doc raw"""
    solveQP_AffineEq

We have a problem of the form:

```math
\underset{x}{\text{minimize}} \quad \frac{1}{2} x^⊤Q x + p^⊤x
```
```math
\text{subject to} \quad Ax = b
```

This yields the KKT system:

```math
[Q A^⊤; A 0] \ [x; λ] = [-p; b]
```

This can be solved with LU Factorization, but we will let Julia's "\" function
handle that.

"""
function solveQP_AffineEq(Q, p, A, b)
    hMat = [Q A'; A spzeros(size(A, 1), size(A, 1))]
    dMat = [-p; b]

    return hMat \ Array(dMat)
end

function solveQP_AffineEq(obj::objectiveQP, constraints::AL_AffineEquality)
    Q = obj.Q
    p = obj.p

    A = constraints.A
    b = constraints.b

    return solveQP_AffineEq(Q, p, A, b)
end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
# wrappers for LLVM-specific functionality

@inline trap() = ccall("llvm.trap", llvmcall, Cvoid, ())
@inline debugtrap() = ccall("llvm.debugtrap", llvmcall, Cvoid, ())

@inline assume(cond::Bool) =
    Base.llvmcall(("declare void @llvm.assume(i1)",
                    "%cond = icmp eq i8 %0, 1
                     call void @llvm.assume(i1 %cond)
                     ret void"), Nothing, Tuple{Bool}, cond)


text:
mutable struct Alignment
    names::Array{String,1} # sequence names
    M::Array{Char,2}
    repr_n:: Array{Int32,1} #number of epresentated species/sequiences - used for
    name_sequence_map:: Dict{String,Int32}

    function Alignment(namesini::Array{String,1}, M::Array{Char,2})
        repr_n = fill(Int32(1), size(M,1))
        names = fill("", size(M,1))
        for i in eachindex(namesini)
            parts = split(namesini[i], ";")
            names[i] = parts[1]
            repr_n[i] = parse(Int32, split(parts[2], "=")[2])
        end
        return Alignment(names, M, repr_n)
    end

    function Alignment(names::Array{String,1}, M::Array{Char,2}, repr_n::Array{Int32,1})
        name_sequence_map = Dict{String,Int32}()
        for i in eachindex(names)
            name_sequence_map[names[i]] = i
        end
        return new(names, M, repr_n, name_sequence_map)
    end

    #primary=true  when the initial alignment is beeing read in (as a reference)
    function  Alignment(file_name::String, allowed_symbols::Set{Char}=Set{Char}(['A','T','G','C','N','-']))::Alignment
        sequences = Array{String,1}()
        names = Array{String,1}()
        num = 0
        lengths = Array{Int32,1}()
        skiped = 0
        println(stderr, "Reading aligment $file_name\n")
        FastaReader(file_name) do fr
            for (n, seq) in fr
                seq = uppercase(seq)
                proper = true
                for s in seq
                    proper = false
                    for sa in allowed_symbols
                        if s == sa
                            proper = true
                            break
                        end
                    end
                    if !proper
                        break
                    end
                end
                if !proper
                    skiped += 1
                    continue
                end
                num += 1
                if mod(num, 10) == 0
                    print(stderr,"Parsed $num sequences \r")
                end
                push!(names, n)
                push!(sequences, seq)
                push!(lengths, length(seq))
            end
        end
        len = maximum(lengths)
        #create alignment array
        m=fill('-', num, len)
        for i in 1:num
            m[i,1:length(sequences[i])] = collect(sequences[i])
        end
        println(stderr)
        println(stderr, "Reading alignment finished, $skiped sequences were skipped due to atypical symbols")
        println(stderr)
        return Alignment(names, m)
    end
end

#sub alignment  - column slice
function sub_alignment(aln::Alignment, interval::UnitRange{Int64})
    names = aln.names
    repr_n = aln.repr_n
    min = aln.M
    mout = Array{Char}(undef, size(min)[1], length(interval))
    for i in 1:size(min)[1]
        mout[i,:] = min[i, interval]
    end
    return(Alignment(names, mout, repr_n))
end

#sub alignment  - subset on particular names
function sub_alignment(aln::Alignment, names::Array{String,1})
    new_names = Array{String,1}() # sequence names
    new_M = Array{Char}(undef, length(names), size(aln.M)[2])
    new_repr_n = Array{Int32,1}() #number of epresentated species/sequiences - used for
    ct = 0
    for n in names
        ct += 1
        ind = aln.name_sequence_map[n]
        push!(new_names, aln.names[ind])
        new_M[ct, :] = aln.M[ind, :]
        push!(new_repr_n, aln.repr_n[ind])
    end
    return(Alignment(new_names, new_M, new_repr_n))
end

function write_to_fasta(aln::Alignment, out_f_name::String)
    f = open(out_f_name,"w")
    for i in 1:size(aln.M, 1)
        n = aln.names[i]
        rn = aln.repr_n[i]
        outn = "$n;size=$rn"
        s = String(aln.M[i,:])
        print(f, ">$outn\n$s\n")
    end
    close(f)
end

function write_to_fasta(aln::Alignment, interval::UnitRange{Int64}, out_f_name::String;remove_gaps=false)
    f = open(out_f_name, "w")
    for i in 1:size(aln.M, 1)
        n = aln.names[i]
        rn = aln.repr_n[i]
        outn = "$n;size=$rn"
        s = String(aln.M[i,interval])
        if remove_gaps
            s=replace(s, "-"=>"N")
        end
        print(f, ">$outn\n$s\n")
    end
    close(f)
end

 """
 Cluster a reagion in an alignemnt and replace custered fragments with a cluster representative
 """
function vsearh_cluster!(aln::Alignment, interval::UnitRange{Int64}, id::Float64)
    #create temorary file in dev/shm
    bindir = Conda.bin_dir(:vsearch)
    tmpdir = mktempdir("/dev/shm"; prefix="vsearch_", cleanup=false)
    in_fasta = "$tmpdir/input.fasta"
    out_fasta = "$tmpdir/centroids.fasta"
    clusters = "$tmpdir/centroids.uc"
    write_to_fasta(aln, interval, in_fasta, remove_gaps=true)
    cluster = run(`$bindir/vsearch --quiet --sizein --cluster_size $in_fasta --uc $clusters --id $id`)
    clustermap = get_clustering(clusters)
    for i in eachindex(aln.names)
        aln.M[i,interval] = aln.M[aln.name_sequence_map[clustermap[aln.names[i]]], interval]
    end
end

function get_max_cl_length(aln::Alignment)
    #create temorary file in dev/shm
    bindir = Conda.bin_dir(:vsearch)
    tmpdir = mktempdir("/dev/shm"; prefix="vsearch_", cleanup=true)
    lenin = size(aln.M)[2]
    target = size(aln.M)[1]
    nmb = 0
    ln = 0
    for i in 10:lenin
        clusters = "$tmpdir/$i.centroids.uc"
        in_fasta = "$tmpdir/$i.input.fasta"
        write_to_fasta(aln, 1:i, in_fasta, remove_gaps=true)
        cluster = run(`$bindir/vsearch --quiet --sizein --cluster_size $in_fasta --uc $clusters --id 0.96`)
        clustermap = get_clustering(clusters)
        ln = i
        nmb = length(unique(values(clustermap)))
        if nmb == target
            break
        end
    end
    println("$nmb serquences at $ln length")
end

function vsearh_cluster!(aln::Alignment, window::Int64, step::Int64, id::Float64)
    lenaln = size(aln.M)[2]
    final = lenaln - window + 1
    starts = collect(1:step:final)
    nb = length(starts)
    ct = 0
    for startsmooth in starts
        ct += 1
        startl = startsmooth
        endl = startsmooth + window -1
        if ct == nb
            endl = lenaln
        end
        vsearh_cluster!(aln, startl:endl, id)
    end
end

"Patses uclust file and collects mapping"
function get_clustering(ucinfile::String)
    out = Dict{String,String}()
    curr = " "
    for l in readlines(ucinfile)
        if l[1] == 'H'
            parts = split(l)
            mem, rep = split(parts[9],";")[1], split(parts[10],";")[1]
            out[mem]  =rep
        end
        if l[1] == 'C'
            parts = split(l)
            mem = split(parts[9], ";")[1]
            out[mem] = mem
        end
    end
    return(out)
end

function char_array_entropy(s::String;debug=false)
    return char_array_entropy(collect(s), debug=debug)
end

function char_array_entropy(s::Array{Char,1}, c::Array{Int,1};debug=false)
    if length(s) != length(c)
        error("Number of symbols and their counts do not match")
    end
    #not perfect way...just generate an array
    chars = Array{Char,1}()
    for i in 1:length(s)
        part=fill(s[i], c[i])
        for cc in part
            push!(chars, cc)
        end
    end
    return(char_array_entropy(chars, debug=debug))
end

function char_array_entropy(s::Array{Char,1};debug=false)
    c = counter(s)
    n = length(s)
    e=0
    consensus = ' '
    ct_gap = 0
    for k in keys(c)
    p = c[k]/n
    e += -1*p*log2(p)
    end
    frac_gap = c['-']/n
    s = (2-e)-(2*frac_gap)
    if debug
        println((2-e), " ", frac_gap)
    end
    consensus = sort(collect(c), by=x->x[2])[1][1]
    return(s, frac_gap,consensus)
end

function aln_conservation(aln::Alignment; ignore_amounts = false)
    entropies = Array{Float64,1}()
    sz = size(aln.M)
    for i in 1:sz[2]
        chars = aln.M[:,i]
        cnts = aln.repr_n
        if ignore_amounts
            entr = char_array_entropy(chars)
        else
            entr = char_array_entropy(chars, Array{Int64,1}(cnts))
        end
        push!(entropies,entr[1])
    end
    return(entropies)
end

"Identfies sequences which are recombinations of others and has no unique snps"
function eval_recomb(aln::Alignment)
    out_dataframe=DataFrame(Read_name = Array{String,1}(),
                Unique_snps = Array{Int64,1}(),
                Maximum_common_fraction = Array{Float64,1}(),
                Maximum_common_fraction_names = Array{String,1}())

    aln_c = aln_conservation(aln)
    m_sub_ini = aln.M[:, aln_c .< 2]
    m_sub = copy(m_sub_ini)
    n_seqs = size(aln.M)[1]
    seqs_length = size(aln.M)[2]
    fractions = Dict{String,Array{Float64,1}}()
    fractions_names = Dict{String,Array{String,1}}()
    for i in 1:n_seqs
        fractions[aln.names[i]] = Array{Float64,1}()
        fractions_names[aln.names[i]] = Array{String,1}()
        for j in 1:n_seqs
            if i != j && (aln.repr_n[i] <  aln.repr_n[j])
                #collect fraction of overlaping snps independant for each reference
                replaced,frac = remove_matching(m_sub_ini[i,:], m_sub_ini[j,:])
                push!(fractions[aln.names[i]], frac)
                push!(fractions_names[aln.names[i]], aln.names[j])
                #collect comparisions
                replaced_acc2,frac2 = remove_matching(m_sub[i,:], m_sub[j,:])
                m_sub[i,:] = replaced_acc2
            end
        end
    end
    out_nonunic = Array{Int64,1}()
    out_max_frac = Array{Float64,1}()
    out_max_frac_names = Array{Float64,1}()
    for i in 1:n_seqs
        nb_non_unic = 0
        if '-' in m_sub[i,:]
            nb_non_unic = counter(m_sub[i,:])['-']
        end
        max_frac_common = 0
        max_frac_common_name = ""
        if length(fractions[aln.names[i]]) > 0
            max_frac_common = maximum(fractions[aln.names[i]])
            max_frac_common_name = join( fractions_names[aln.names[i]][ fractions[aln.names[i]] .== max_frac_common ], " ")
        end
        nb_unic = length(m_sub[i,:]) - nb_non_unic
        push!(out_dataframe, [aln.names[i], nb_unic, max_frac_common, max_frac_common_name])
    end
    return(out_dataframe)
end

"changes a symbol in act array to '' if in the ref array is the same symbol"
function remove_matching(act::Array{Char,1}, ref::Array{Char,1})
    out = copy(act)
    ct = 0
    for i in 1:length(act)
        if act[i] == ref[i]
            out[i] = '-'
            ct += 1
        end
    end

    return(out, ct/length(act))
end


text:
module StructuredMatrices

using SIMDPirates, SLEEFPirates, VectorizationBase, LinearAlgebra, PaddedMatrices, LoopVectorization, StackPointers
using PaddedMatrices: AbstractFixedSizeVector, AbstractFixedSizeMatrix
using VectorizationBase: REGISTER_SIZE, REGISTER_COUNT
using ReverseDiffExpressionsBase
using Parameters

import PaddedMatrices: param_type_length, type_length

export addmul!, submul!, inv′, ∂inv′,
        UpperTriangularMatrix, LowerTriangularMatrix, SymmetricMatrixL,
        AutoregressiveMatrixLowerCholeskyInverse, AutoregressiveMatrix,
        MutableLowerTriangularMatrix, MutableUpperTriangularMatrix

@noinline ThrowBoundsError(args...) = throw(BoundsError(args...))
@inline binomial2(n::Integer) = (n*(n-1)) >>> 1

# include("static_ranges.jl")
include("triangular_representation_utilities.jl")
include("symmetric_matrix.jl")
include("triangular_matrix.jl")
include("triangle_inverse.jl")
include("vector_of_triangular_matrix_operations.jl")
include("autoregressive_matrix.jl")
include("block_diagonal.jl")
include("adjoint.jl")
include("decompositions.jl")
include("triangular_equations.jl")
include("ragged_matrix.jl")

@def_stackpointer_fallback rank_update rank_update! reverse_cholesky_grad ∂rank_update
include("precompile.jl")
_precompile_()
function __init__()
    @add_stackpointer_method rank_update rank_update! reverse_cholesky_grad ∂rank_update
    _precompile_()
end

end # module


text:
# DESCRIPTION
# Model order reduction based on residualization.
#
# INPUT
# - A: system matrix
# - B: control matrix
# - idxss: indices for slow states
# - idxfs: indices fro fast states
#
# OUTPUT
# - Ared: reduced-order system matrix
# - Bred: reduced-order control matrix

function ModRed(A,B,idxss,idxfs)
    Ared=A[idxss,idxss]-A[idxss,idxfs]/A[idxfs,idxfs]*A[idxfs,idxss]
    Bred=B[idxss,1:4]-A[idxss,idxfs]/A[idxfs,idxfs]*B[idxfs,1:4]
    return Ared, Bred
end

max_stars_repo_name: umbertosaetti/J-GenHel

max_issues_repo_name: umbertosaetti/J-GenHel

max_forks_repo_name: umbertosaetti/J-GenHel


text:
using Documenter
using Rifraf

makedocs(
    modules = [Rifraf],
    doctest = false
)


text:
# Scalar operations

# @inline Base.@pure vshr(a::Int64, b) = llvmcall("%res = ashr i64 %0, %1\nret i64 %res", Int64, Tuple{Int64,Int64}, a, b % Int64)
# @inline Base.@pure vshr(a::Int32, b) = llvmcall("%res = ashr i32 %0, %1\nret i32 %res", Int32, Tuple{Int32,Int32}, a, b % Int32)

# @inline Base.@pure vshr(a::Int16, b) = llvmcall("%res = ashr i16 %0, %1\nret i16 %res", Int16, Tuple{Int16,Int16}, a, b % Int16)
# @inline Base.@pure vshr(a::Int8, b) = llvmcall("%res = ashr i8 %0, %1\nret i8 %res", Int8, Tuple{Int8,Int8}, a, b % Int8)

# @inline Base.@pure vshr(a::UInt64, b) = llvmcall("%res = lshr i64 %0, %1\nret i64 %res", UInt64, Tuple{UInt64,UInt64}, a, b % UInt64)
# @inline Base.@pure vshr(a::UInt32, b) = llvmcall("%res = lshr i32 %0, %1\nret i32 %res", UInt32, Tuple{UInt32,UInt32}, a, b % UInt32)

# @inline Base.@pure vshr(a::UInt16, b) = llvmcall("%res = lshr i16 %0, %1\nret i16 %res", UInt16, Tuple{UInt16,UInt16}, a, b % UInt16)
# @inline Base.@pure vshr(a::UInt8, b) = llvmcall("%res = lshr i8 %0, %1\nret i8 %res", UInt8, Tuple{UInt8,UInt8}, a, b % UInt8)

function binary_op(op, W, @nospecialize(UNUSED::Type{T}), W2 = W) where {T}
    @assert W == W2
    ty = LLVM_TYPES[T]
    if isone(W)
        V = T
    else
        ty = "<$W x $ty>"
        V = NTuple{W,VecElement{T}}
    end
    instrs = "%res = $op $ty %0, %1\nret $ty %res"
    call = :(llvmcall($instrs, $V, Tuple{$V,$V}, data(v1), data(v2)))
    W > 1 && (call = Expr(:call, :Vec, call))
    Expr(:block, Expr(:meta, :inline), call)
end
# @generated function binary_operation(::Val{op}, v1::V1, v2::V2) where {op, V1, V2}
#     M1, N1, W1, T1 = description(V1)
#     M2, N2, W2, T2 = description(V2)

#     lc = Expr(:call, :llvmcall, join(instrs, "\n"), )
# end
# function binary_op(op, W, @nospecialize(_::Type{T})) where {T}
#     ty = 'i' * string(8*sizeof(T))
#     binary_op(op, W, T, ty)
# end

# Integer
    # vop = Symbol('v', op)
for (op,f) ∈ [("add",:+),("sub",:-),("mul",:*),("shl",:<<)]
    ff = Symbol('v', op); _ff = Symbol(:_, ff)
    @eval @generated $_ff(v1::T, v2::T) where {T<:Integer} = binary_op($op * (T <: Signed ? " nsw" : " nuw"), 1, T)
    @eval @generated $_ff(v1::Vec{W1,T}, v2::Vec{W2,T}) where {W1,W2,T<:Integer} = binary_op($op * (T <: Signed ? " nsw" : " nuw"), W1, T, W2)
    @eval @generated Base.$f(v1::Vec{W1,T}, v2::Vec{W2,T}) where {W1,W2,T<:Integer} = binary_op($op, W1, T, W2)
    @eval Base.@pure @inline $ff(v1::T, v2::T) where {T} = $_ff(v1, v2)
    @eval @inline $ff(v1, v2) = ((v3, v4) = promote(v1, v2); $ff(v3, v4))
end
for (op,f) ∈ [("div",:÷),("rem",:%)]
    ff = Symbol('v', op); _ff = Symbol(:_, ff)
    @eval @generated Base.$f(v1::Vec{W1,T}, v2::Vec{W2,T}) where {W1,W2,T<:Integer} = binary_op((T <: Signed ? 's' : 'u') * $op, W1, T, W2)
    @eval @generated $_ff(v1::T, v2::T) where {T<:Integer} = binary_op((T <: Signed ? 's' : 'u') * $op, 1, T)
    @eval @inline $_ff(v1::AbstractSIMD, v2::AbstractSIMD) = Base.$f(v1, v2)
    @eval Base.@pure @inline $ff(v1::T, v2::T) where {T} = $_ff(v1, v2)
    @eval @inline $ff(v1, v2) = ((v3, v4) = promote(v1, v2); $ff(v3, v4))
end
@inline vcld(x, y) = vadd(vdiv(vsub(x,one(x)), y), one(x))
@inline function vdivrem(x, y)
    d = vdiv(x, y)
    r = vsub(x, vmul(d, y))
    d, r
end
for (op,f,s) ∈ [("ashr",:>>,0x01),("lshr",:>>,0x02),("lshr",:>>>,0x03),("and",:&,0x03),("or",:|,0x03),("xor",:⊻,0x03)]
    _ff = Symbol(:_, 'v', op)
    fdef = Expr(:where, :(Base.$f(v1::Vec{W1,T}, v2::Vec{W2,T})), :W1, :W2)
    ffdef = Expr(:where, :($_ff(v1::T, v2::T)))
    if iszero(s & 0x01)
        push!(fdef.args, :(T <: Unsigned))
        push!(ffdef.args, :(T <: Unsigned))
    elseif iszero(s & 0x02)
        push!(fdef.args, :(T <: Signed))
        push!(ffdef.args, :(T <: Signed))
    else
        push!(fdef.args, :T)
        push!(ffdef.args, :T)
    end
    @eval @generated $fdef = binary_op($op, W1, T, W2)
    @eval @generated $ffdef = binary_op($op, 1, T)
    if op === "ashr" # replace first iteration with <<
        f = :(<<)
    end
    @eval begin
        @inline Base.$f(v::Vec, i::Integer) = ((x, y) = promote(v, i); $f(x, y))
        @inline Base.$f(i::Integer, v::Vec) = ((x, y) = promote(i, v); $f(x, y))
    end
end
for (op,f) ∈ [("lshr",:>>),("ashr",:>>),("and",:&),("or",:|),("xor",:⊻)]
    ff = Symbol('v', op); _ff = Symbol(:_, ff)
    @eval Base.@pure @inline $ff(v1::T, v2::T) where {T} = $_ff(v1, v2)
    @eval @inline $ff(v1, v2) = ((v3, v4) = promote(v1, v2); $ff(v3, v4))
end

for (op,f,ff) ∈ [("fadd",:+,:vadd),("fsub",:-,:vsub),("fmul",:*,:vmul),("fdiv",:/,:vdiv),("frem",:%,:vrem)]
    @eval begin
        @generated Base.$f(v1::Vec{W1,T}, v2::Vec{W2,T}) where {W1,W2,T<:Union{Float32,Float64}} = binary_op($(op * " nsz arcp contract afn reassoc"), W1, T, W2)
        @generated $ff(v1::Vec{W1,T}, v2::Vec{W2,T}) where {W1,W2,T<:Union{Float32,Float64}} = binary_op($op, W1, T, W2)
        @inline $ff(a::T, b::T) where {T <: Union{Float32,Float64}} = $f(a, b)
    end
end

@inline Base.:(/)(a::Vec{W,<:Integer}, b::Vec{W,<:Integer}) where {W} = float(a) / float(b)

for op ∈ [:⊻, :&, :|]
    @eval begin
        @inline function Base.$op(v1::AbstractSIMDVector{W}, v2::AbstractSIMDVector{W}) where {W}
            v3, v4 = promote(v1, v2)
            $op(v3, v4)
        end
    end
end

function promote_shift_quote(op::Symbol, ::Type{T1}, ::Type{T2}) where {T1, T2}
    s1 = sizeof(T1); s2 = sizeof(T2);
    @assert s1 != s2 "op: $op, T1 === $T1, T2 === $T2"
    if s1 < s2
        newT = T1 <: Signed ? signed(T2) : unsigned(T2)
        f = Expr(:call, op, Expr(:call, :convert, newT, :v1), :v2)
    else
        newT = T2 <: Signed ? signed(T1) : unsigned(T1)
        f = Expr(:call, op, :v1, Expr(:call, :convert, newT, :v2))
    end
    Expr(:block, Expr(:meta, :inline), f)
end
@inline function Base.:(<<)(v1::AbstractSIMDVector{W,T1}, v2::AbstractSIMDVector{W,T2}) where {W,T1<:SignedHW,T2<:UnsignedHW}
    convert(T1, convert(T2, v1) << v2)
end
@inline function Base.:(<<)(v1::AbstractSIMDVector{W,T1}, v2::AbstractSIMDVector{W,T2}) where {W,T1<:UnsignedHW,T2<:SignedHW}
    convert(T1, convert(T2, v1) << v2)
end
@inline function Base.:(>>)(v1::AbstractSIMDVector{W,T1}, v2::AbstractSIMDVector{W,T2}) where {W,T1<:SignedHW,T2<:UnsignedHW}
    convert(T1, convert(T2, v1) >> v2)
end
@inline function Base.:(>>)(v1::AbstractSIMDVector{W,T1}, v2::AbstractSIMDVector{W,T2}) where {W,T1<:UnsignedHW,T2<:SignedHW}
    convert(T1, convert(T2, v1) >> v2)
end
@inline function Base.:(>>>)(v1::AbstractSIMDVector{W,T1}, v2::AbstractSIMDVector{W,T2}) where {W,T1<:SignedHW,T2<:UnsignedHW}
    convert(T1, convert(T2, v1) >>> v2)
end
@inline function Base.:(>>>)(v1::AbstractSIMDVector{W,T1}, v2::AbstractSIMDVector{W,T2}) where {W,T1<:UnsignedHW,T2<:SignedHW}
    convert(T2, v1 >>> convert(T1, v2))
end
for op ∈ [:(<<), :(>>), :(>>>), :(&), :(|)]
    @eval begin
        @generated function Base.$op(v1::AbstractSIMDVector{W,T1}, v2::AbstractSIMDVector{W,T2}) where {W,T1<:IntegerTypes,T2<:IntegerTypes}
            promote_shift_quote($(QuoteNode(op)), T1, T2)
        end
    end
end


text:
ENV["GKSwstype"] = "nul"  # For StatsPlots

import Pkg; Pkg.activate("../../../")
include("../../info.jl")
include("scenarios.jl")

using CytofDiffDensity; const cdd = CytofDiffDensity
using Util
import Random
using DrWatson
using MCMC
using BSON
using LaTeXStrings
using Distributions  # required for `BSON.load`
using StatsPlots
using StatsFuns

# Plot settings
plotsize = (450, 450)
Plots.scalefontsizes()
Plots.scalefontsizes(1.5)

simname = "v3"
make_resultsdir(sim) = mkpath("$(Info.resultsdir_simstudy)/$(simname)/$(savename(sim))")

function plot_simdata!(simdata; ygrid=nothing)
  histogram!(simdata.yC, normalize=true, label=nothing, la=0, alpha=.3, color=:blue)
  histogram!(simdata.yT, normalize=true, label=nothing, la=0, alpha=.3, color=:red)
  ygrid === nothing && (ygrid = cdd.make_ygrid([simdata.yC; simdata.yT], 100))
  plot!(ygrid, pdf.(simdata.mmC, ygrid), label=nothing, color=:blue)
  plot!(ygrid, pdf.(simdata.mmT, ygrid), label=nothing, color=:red)
end

function run(sim)
  sim = (; sim...)

  # Simulate data.
  simdata = scenarios(sim.snum, Ni=100_000, seed=1)

  # Make results dir if needed.
  resultsdir = make_resultsdir(sim)

  # Create Priors.
  y = [simdata.yC; simdata.yT]

  # Plot and save data
  imgdir = mkpath(joinpath(resultsdir, "img"))
  plot(size=plotsize)
  plot_simdata!(simdata)
  savefig(joinpath(imgdir, "data.pdf"))
  closeall()

  Random.seed!(0)
  mu_prior = Normal(mean(y), std(y))  # independent priors for mu_k.
  if sim.skewtmix
    # SkewT mixture
    model = Gtilde(simdata.yC, simdata.yT, sim.K, true, mu=mu_prior, skew=true, tdist=true)
  else
    # Normal mixture
    model = Gtilde(simdata.yC, simdata.yT, sim.K, true, mu=mu_prior, skew=false, tdist=false)
  end
  exclude = [:vC, :vT, :zetaC, :zetaT, :lambdaC, :lambdaT]

  init = MCMC.make_init_state(model)
  spl = make_sampler(model, init=init)

  nburn, nsamps, thin = (20000, 4000, 2)
  callback = make_callback(model, nburn=nburn, nsamps=nsamps, thin=thin)

  Util.redirect_stdout_to_file(joinpath(resultsdir, "log.txt")) do
    # Run MCMC.
    chain, metrics = mcmc(spl, nsamps, init=init, nburn=nburn, thin=thin,
                          callback=callback, exclude=exclude);

    # Save chain and metrics.
    BSON.bson(joinpath(resultsdir, "results.bson"),
              Dict(:chain => chain, :metrics => metrics, :simdata => simdata,
                   :model => model))
  end
end

function _postprocess(sim, resultsdir)
  sim = (; sim...)

  # Make image dir if needed.
  imgdir = mkpath(joinpath(resultsdir, "img"))

  # Load results.
  # chain, metrics, simdata, model.
  r = (; BSON.load(joinpath(resultsdir, "results.bson"))...)

  # Print summary stats
  cdd.printsummary(r.chain, r.metrics); println()

  # Print model info
  cdd.print_model_info(r.model)

  # Plot loglike
  plot(r.metrics[:loglike], label=nothing)
  savefig(joinpath(imgdir, "loglike.pdf"))
  closeall()

  # Plot posterior for each multivariate parameter.
  for sym in [:mu, :sigma, :nu, :phi, :etaC, :etaT]
    p = hcat(getindex.(r.chain, sym)...)'

    plot(p, label=nothing)
    savefig(joinpath(imgdir, "$(sym)-trace.pdf"))
    closeall()

    boxplot(p, label=nothing)
    savefig(joinpath(imgdir, "$(sym).pdf"))
    closeall()
  end

  # Plot posterior for each univariate parameter.
  for sym in [:tau]
    p = vec(hcat(getindex.(r.chain, sym)...))

    plot(p, label=nothing)
    savefig(joinpath(imgdir, "$(sym)-trace.pdf"))
    closeall()

    histogram(p, label=nothing, normalize=true)
    savefig(joinpath(imgdir, "$(sym).pdf"))
    closeall()
  end

  # Plot posterior density for C and T.
  y = [r.simdata.yC; r.simdata.yT]
  ygrid = cdd.make_ygrid(y, 100)
  plot(size=plotsize)
  cdd.plot_post_density!(r.chain, ygrid)
  cdd.plot_simtruth(r.simdata.mmC, r.simdata.mmT, ygrid, label=nothing)
  savefig(joinpath(imgdir, "post-density.pdf"))
  closeall()

  # Plot gammas (zero part)
  begin
    pzero = Pzero(NC=r.simdata.NC, NT=r.simdata.NT,
                  QC=r.simdata.QC, QT=r.simdata.QT)
    pzero_inference = cdd.infer_Pzero1(pzero, 10000)
    gammaC_post = pzero_inference.distC
    gammaT_post = pzero_inference.distT
    gamma_grid = cdd.make_gamma_grid(gammaC_post, gammaT_post, a=.001)

    plot(size=plotsize)
    plot!(gamma_grid, pdf.(gammaC_post, gamma_grid), label=nothing, color=:blue)
    plot!(gamma_grid, pdf.(gammaT_post, gamma_grid), label=nothing, color=:red)
    cdd.plot_gamma_uq!(gammaC_post, color=:blue, Q=pzero.QC, N=pzero.NC,
                       truth=r.simdata.gammaC)
    cdd.plot_gamma_uq!(gammaT_post, color=:red, Q=pzero.QT, N=pzero.NT,
                       truth=r.simdata.gammaT)
    savefig(joinpath(imgdir, "gamma.pdf"))
    closeall()
  end

  # Plot F̃ᵢ
  begin
    # Plot Fi-tilde
    plot(size=plotsize)
    cdd.plot_Fi_tilde_cdf!(r.model, r.chain, pzero, N=2000)
    ylims!(0, 1)
    xlabel!(L"\tilde{y}")
    ylabel!("CDF")
    savefig(joinpath(imgdir, "Fi-tilde-postmean.pdf"))
    closeall()

    # Plot Fi
    plot(size=plotsize)
    cdd.plot_Fi_tilde_cdf!(r.model, r.chain, pzero, exponentiate=true, N=2000)
    ylims!(0, 1)
    xlabel!("y")
    ylabel!("CDF")
    savefig(joinpath(imgdir, "Fi-postmean.pdf"))
    closeall()
  end
end

function postprocess(sim)
  # Make results dir if needed.
  resultsdir = make_resultsdir(sim)

  Util.redirect_stdout_to_file(joinpath(resultsdir, "info.txt")) do
    _postprocess(sim, resultsdir)
  end
end

function parse_dic(model_info::String)
  m = match(r"(?<=DIC:\s)\d+\.\d+", model_info).match
  return parse(Float64, m)
end