text:
require("Feature.jl")
require("Eisner.jl")

immutable Model
	weights::Vector{Real}
	featdict::Dict{String, Int}
end

immutable Word
	id::Int
	form::String
	lemma::String
	cpostag::String
	postag::String
	feats::Array
	head::Int
	deprel::String
end

function readconlldata(fpath)
	info("Reading $fpath ...")
	const file = open(fpath, "r")

	sentences = Array{Word,1}[] # a list of sent
	sent = Word[] # a list of Word
	for line in eachline(file)
		if line != "\n"
			cols = split(line, "\t")
			token = Word(
				int(cols[1]), 			# id
				cols[2], 				# form
				cols[3], 				# lemma
				cols[4], 				# cpostag
				cols[5], 				# postag
				split(cols[6],"|"), 	# feats
				int(cols[7]), 			# head
				cols[8] 				# deprel
			)
			push!(sent, token)
		else
			push!(sentences, sent)
			sent = Word[]
		end
	end
	close(file)

	info("Total number of sentences: $(length(sentences))")
	sentences
end

function getedgescores(edge_feats, weights::Vector{Real}, n::Int)
	# n: number of words in sentence
	# n = sqrt(length(edge_feats)) # total number of edges = n*n
	edge_scores = Dict{(Int, Int), Real}()
	sizehint(edge_scores, n*n) # n*n: number of edges to consider

	for index_head = 0:n
		for index_mod = 1:n
			if index_head != index_mod
				score = 0
				for findex in edge_feats[(index_head, index_mod)]
					score += weights[findex]
				end
			end
			edge_scores[(index_head, index_mod)] = score
		end
	end

	edge_scores
end


text:
# https://github.com/JuliaMusic/MusicXML.jl/issues/46
using MusicXML

@importMX
xml_note = parsemusicxml_partial("""
<note>
  <grace/>
  <pitch>
    <step>G</step>
    <octave>4</octave>
  </pitch>
  <voice>1</voice>
  <type>16th</type>
  <stem>up</stem>
  <beam number="1">end</beam>
  <beam number="2">end</beam>
</note>
""")

Note(xml_note)

ext: jl


text:
using LinearAlgebra
using FermiCG
using Printf
using Test
using LinearMaps
using Arpack
using Random
using Profile

@testset "symm" begin

    N = 1000
    Nv = N*(N+1)÷2
    println(N)
    println(Nv)
    A = Diagonal(rand(N)) .- .5 + 1.0*rand(N,N)
    A = A'+A
    Hv = FermiCG.SymDenseMat{Float64}(A)

    Hm = Matrix(Hv)
    v = rand(Float64,N)

    r1 = A*v
    r2 = Hv*v
    r3 = Hm*v

    @test all(isapprox(r1,r2, atol=1e-10))
    @test all(isapprox(r1,r3, atol=1e-10))

    nr = 2

    e1,v1 = Arpack.eigs(A, nev = nr, which=:SR)
    e2,v2 = Arpack.eigs(Hv, nev = nr, which=:SR)

    @test all(isapprox(e1,e2, atol=1e-10))
end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
using ForecastData
using Test

using DataFrames

@testset "ForecastData.jl" begin

    # co2
    data = co2()
    @test data isa DataFrame
    @test size(data) == (4609,2)

    data = co2("full")
    @test data isa DataFrame
    @test size(data) == (17166,17)

    # london
    data = london()
    @test data isa DataFrame
    @test size(data) == (132,15)

    # quakes
    data = quakes()
    @test data isa DataFrame
    @test size(data) == (181,2)

    data = quakes("full")
    @test data isa DataFrame
    @test size(data) == (11924,22)

    #seaborne
    data = seaborne()
    @test data isa DataFrame
    @test size(data) == (2199,4)

    data = seaborne("full")
    @test data isa DataFrame
    @test size(data) == (13209,10)

end


text:
using Pkg; Pkg.activate("."); Pkg.instantiate(); #install required packages
server = "ai.ku.edu.tr/"
if length(ARGS)==0
    println("clevr home folder is not specified. Pre-processed data will be downloaded from the servers(70GB) if they are not in the data folder? (yes or no)")
    if readline() == "yes"
        !isfile("data/train.bin")  && download(server*"data/mac-network/train.bin","data/train.bin")
        !isfile("data/val.bin")    && download(server*"data/mac-network/val.bin","data/val.bin")
        !isfile("data/train.json") && download(server*"data/mac-network/train.json","data/train.json")
        !isfile("data/val.json")   && download(server*"data/mac-network/val.json","data/val.json")
        !isfile("data/dic.json")   && download(server*"data/mac-network/dic.json","data/dic.json")
        exit()
    else
        error("No data available!")
    end
else
    CLEVR_HOME = ARGS[1]
    println(CLEVR_HOME)
end

include("preprocess.jl")
println("Starting to preprocess CLEVR questions @: $(CLEVR_HOME)/questions ")
println("Training and validation questions will be processed...")
if isfile("data/train.json") && isfile("data/val.json")
    println("You already had processed questions in your data folder. Do you want to re-process the data? (yes or no)")
    readline() == "yes" && preprocess([CLEVR_HOME])
else
    preprocess([CLEVR_HOME])
end

include("extract_features.jl")
println("Starting to preprocess CLEVR images @: $(CLEVR_HOME)/images ")
println("Training and validation images will be processed...")

if isfile("data/train.bin") && isfile("data/val.bin")
    println("You already had processed images in your data folder. Do you want to re-process the data? (yes or no)")
    readline() == "yes" && extract_features([CLEVR_HOME])
else
    extract_features([CLEVR_HOME])
end

max_stars_repo_name: feyzaakyurek/MAC-Network

max_issues_repo_name: feyzaakyurek/MAC-Network

max_forks_repo_name: feyzaakyurek/MAC-Network


text:
#
# Tree hierarchy
#

using JSON

abstract NodeType
abstract DirectoryNode <: NodeType
abstract FileNode <: NodeType
abstract FunctionNode <: NodeType

const typenames = Dict{Type,Symbol}(
    DirectoryNode   => :directory,
    FileNode        => :file,
    FunctionNode    => :function
)

type Node{T<:NodeType}
    props::Dict{Symbol,Any}
    nodes::Dict{String,Node}

    Node() = new(Dict{Symbol,Any}(), Dict{String,Node}())
end

function JSON.lower{T<:NodeType}(node::Node{T})
    return Dict{Symbol, Any}(
        :type => typenames[T],
        :nodes => node.nodes,
        node.props...
    )
end

Base.getindex(node::Node, key::Symbol) = node.props[key]
Base.setindex!(node::Node, val, key::Symbol) = (node.props[key] = val)

#
# Helper methods
#

function create_path!(root::Node{DirectoryNode}, path::String)
    dirpath, filename = splitdir(path)
    dirs = split(dirpath, '/')
    parent = root
    for dir in dirs
        parent = get!(parent.nodes, dir, Node{DirectoryNode}())
    end

    return parent, filename
end

max_stars_repo_name: OhioAdam/Data-Tree-Modeler

max_issues_repo_name: OhioAdam/Data-Tree-Modeler

max_forks_repo_name: OhioAdam/Data-Tree-Modeler


text:
# This file defines the Condition type, used to define the extra conditions on
#   transitions.

# It is important to note that conditions that define an inequality relation,
#   the relation is ALWAYS stated as
#     "attribute operator value".
# For example, if something requires a physical fitness score of 14.0, the
#   condition is written as "fitness" >= 14.0

# The Condition type requires no extra types.
requiredTypes = [ ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( typePath, reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end  # for reqType in requiredTypes

# This type is not exported.
"""
This type defines a condition on a single attribute of a personnel member, to be
used for checking if a transition will take place.

The type contains the following fields:
* `attr::String`: the attribute upon which the prerequisite is defined.
* `val::T <: Union{Real, String, Vector{String}}`: the value with which the
  attribute is compared.
* `rel::Function`: the relation between the attribute and the value. Valid
relations are `==`, `!=`, `∈`, and `∉`.

A condition is always of the form `attr rel val`.
"""
type Condition{T <: Union{Real, String, Vector{String}}}

    attr::String
    val::T
    rel::Function

    # Constructor.
    function Condition{T}( attr::String, rel::Function, val::T ) where T

        if ( T <: Real ) && ( rel ∉ [ ==, !=, <, <=, >, >= ] )
            error( "A numeric value permits only ==, !=, <, <=, >, and >= as operators." )
        elseif ( T === String ) && ( rel ∉ [ ==, != ] )
            error( "A string value permits only == and != as operators." )
        elseif ( T === Vector{String} ) && ( rel ∉ [ ∈, ∉ ] )
            error( "A vector of string values permits only ∈ and ∉ as operators." )
        end  # if ( T <: Real ) && ...

        newPrereq = new()
        newPrereq.attr = attr
        newPrereq.val = val
        newPrereq.rel = rel
        return newPrereq

    end  # Condition( attr, rel, val )

end  # type Condition


text:
# Script for creating Figure 6, 7, 8, 9 in scale manuscript

using IcanProj
using DataFrames
using JFSM2
using PyPlot
using PyCall
using Statistics

function compute_error(path, res_coarse, variable)

    experiment = 1:32

    res = Vector{Dict}(undef, length(experiment))

    for iexp in experiment

        @show iexp

        # File names

        file_fine = get_filename(path, variable, "1km", iexp)

        file_coarse = get_filename(path, variable, res_coarse, iexp)

        # Load data

        df_links = link_results(file_fine, file_coarse)

        data_coarse, data_aggregated, ngrids = unify_results(file_fine, file_coarse, df_links, variable)

        # Compute metrics per grid cell

        err = data_coarse - data_aggregated

        rmse = sqrt.(mean(err.^2 , dims = 1))

        meanref = mean(data_aggregated, dims = 1)

        meancmp = mean(data_coarse, dims = 1)

        nrmse = rmse ./ meanref

        bias = meancmp .- meanref

        abs_bias = abs.(bias)

        perc_bias = 100*(meancmp ./ meanref .- 1)

        # Return as averages

        weights = ngrids[:]/sum(ngrids[:])

        res[iexp] = Dict("rmse" => sum(rmse[:] .* weights),
                         "meanref" => sum(meanref[:] .* weights),
                         "meancmp" => sum(meancmp[:] .* weights),
                         "nrmse" => sum(nrmse[:] .* weights),
                         "bias" => sum(bias[:] .* weights),
                         "abs_bias" => sum(abs_bias[:] .* weights),
                         "perc_bias" => sum(perc_bias[:] .* weights))

    end

    return res

end

function compute_error_all(path, variable)

    res_coarse = ["5km", "10km", "25km", "50km"]

    resall = Array{Array{Dict,1}}(undef, length(res_coarse))

    for i in 1:length(res_coarse)

        @show res_coarse[i]

        resall[i] = compute_error(path, res_coarse[i], variable)

    end

    return resall

end

function plot_error_scales(resall, variable, ylabel_left, ylabel_right)

data_rmse = zeros(5, 32)
data_bias = zeros(5, 32)

for i in 1:4, j in 1:32
    data_rmse[i+1, j] = resall[i][j]["rmse"]
    data_bias[i+1, j] = resall[i][j]["abs_bias"]
end

df_cfg = cfg_table()

exchng_0 = convert(Array{Bool}, df_cfg[:exchng] .== 0)
exchng_1 = convert(Array{Bool}, df_cfg[:exchng] .== 1)

resolutions = [1, 5, 10, 25, 50]

py"""
import matplotlib.pyplot as plt

fig, ax = plt.subplots(1, 2, figsize=(7, 3.5), tight_layout = True)

ax[0].plot($resolutions, $(data_rmse[:, exchng_0]), color = 'red', label = 'Option 0')
ax[0].plot($resolutions, $(data_rmse[:, exchng_1]), color = 'blue', label = 'Option 1', linestyle = '--')

ax[0].set_ylabel($ylabel_left)

ax[0].set_xticks($resolutions)
ax[0].set_xlabel("Resolution (km)")

ax[0].annotate("(A)", xy=[0.05, 0.9], xycoords='axes fraction', fontsize=12)

ax[1].plot($resolutions, $(data_bias[:, exchng_0]), color = 'red', label = 'Option 0')
ax[1].plot($resolutions, $(data_bias[:, exchng_1]), color = 'blue', label = 'Option 1', linestyle = '--')

ax[1].yaxis.tick_right()

ax[1].yaxis.set_label_position('right')

ax[1].set_ylabel($ylabel_right)

ax[1].set_xticks($resolutions)
ax[1].set_xlabel("Resolution (km)")

ax[1].annotate("(B)", xy=[0.05, 0.9], xycoords='axes fraction', fontsize=12)

plt.show()

plt.savefig($(joinpath(pathfig, "$(variable).pdf")), dpi = 600)
"""

    return nothing

end

# Error different scales - overall settings

pathres = "/data04/jmg/fsm_simulations/netcdf/fsmres_forest"

pathfig = joinpath(dirname(pathof(IcanProj)), "..", "plots", "error_scales")

# Snow water equivalent

variable = "swe"

resall = compute_error_all(pathres, variable)

plot_error_scales(resall, variable, "RMSE (\$kg/m^2\$)", "MAB (\$kg/m^2\$)")

# Net radiation

variable = "rnet"

resall = compute_error_all(pathres, variable)

plot_error_scales(resall, variable, "RMSE (\$W/m^2\$)", "MAB (\$W/m^2\$)")

# Sensible heat fluxes

variable = "hatmo"

resall = compute_error_all(pathres, variable)

plot_error_scales(resall, variable, "RMSE (\$W/m^2\$)", "MAB (\$W/m^2\$)")

# Latent heat fluxes

variable = "latmo"

resall = compute_error_all(pathres, variable)

plot_error_scales(resall, variable, "RMSE (\$W/m^2\$)", "MAB (\$W/m^2\$)")


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# =================================================================================================
"""
    pix_read_tiff(
        filename::AbstractString;
        page::Integer = Int32(1)
    )::Union{Pix, Nothing}

Read a TIFF image from the specified file.  Returns `nothing` on error.

__Arguments:__

| T | Name     | Default    | Description
|:--| :------- | :--------- | :----------
| R | filename |            | The name of the TIFF file to load.
| O | page     | `Int32(1)` | The image to load from the file.

__Restrictions:__

  * `page` - Must be greater than 0.

__Details:__

TIFF files can contain multiple images.  The `page` parameter allows you to specify which image
you want to load.
"""
function pix_read_tiff(
            filename::AbstractString;
            page::Integer = Int32(1)
        )::Union{Pix, Nothing}
    local retval = nothing

    local result = @threadcall(
        (:pixReadTiff, LEPTONICA),
        Ptr{Cvoid},
        (Cstring, Cint),
        filename,
        page - 1
    )

    if result != C_NULL
        retval = Pix(result)
    end

    return retval
end

# =================================================================================================
"""
    pix_read_tiff(
        stream::IO;
        page::Integer = Int32(1)
    )::Union{Pix, Nothing}

Read a TIFF image from the specified stream.  Returns `nothing` on error.

__Arguments:__

| T | Name     | Default    | Description
|:--| :------- | :--------- | :----------
| R | stream   |            | The IO stream to read the TIFF file from.
| O | page     | `Int32(1)` | The image to load from the file.

__Restrictions:__

  * `page` - Must be greater than 0.

__Details:__

This implementation mirrors the API provided by Leptonica when you pass in a FILE pointer. This
assumes that the remainder of the stream contains a TIFF image.

TIFF files can contain multiple images.  The `page` parameter allows you to specify which image
you want to load.
"""
function pix_read_tiff(
            stream::IO;
            page::Integer = Int32(1)
        )::Union{Pix, Nothing}
    local data = read(stream)
    return pix_read_tiff(data; page = page)
end

# =================================================================================================
"""
    pix_read_tiff(
        data::AbstractArray{UInt8};
        page::Integer = Int32(1)
    )::Union{Pix, Nothing}

Read a TIFF image from the byte array.  Returns `nothing` on error.

__Arguments:__

| T | Name | Default    | Description
|:--| :--- | :--------- | :----------
| R | data |            | The byte array to read the TIFF image from.
| O | page | `Int32(1)` | The image to load from the file.

__Restrictions:__

  * `page` - Must be greater than 0.

__Details:__

TIFF files can contain multiple images.  The `page` parameter allows you to specify which image
you want to load.
"""
function pix_read_tiff(
            data::AbstractArray{UInt8};
            page::Integer = Int32(1)
        )::Union{Pix, Nothing}
    local retval = nothing

    local result = ccall(
        (:pixReadMemTiff, LEPTONICA),
        Ptr{Cvoid},
        (Ptr{UInt8}, Csize_t, Cint),
        data,
        length(data),
        page - 1
    )

    if result != C_NULL
        retval = Pix(result)
    end

    return retval
end

# =================================================================================================
"""
    valid_tiff_compression(
        pix::Pix,
        compression::IFF
    )::Bool

Verify that the compression can handle the bit depth of the image.

__Arguments:__

| T | Name        | Default    | Description
|:--| :---------- | :--------- | :----------
| R | pix         |            | The image to check on.
| R | compression |            | The compression being requested by th client.

__Details:__

This method assumes the compression is a valid TIFF compression format.  Some of the TIFF formats
only handle 1bpp images.  If the image is 1bpp then it will return true.  Otherwise it will return
true if the compression supports non 1bpp images.
"""

function valid_tiff_compression(
            pix::Pix,
            compression::IFF
        )::Bool
    if compression == IFF_TIFF
        true
    elseif compression == IFF_TIFF_LZW
        true
    elseif compression == IFF_TIFF_ZIP
        true
    elseif compression == IFF_TIFF_JPEG
        true
    elseif  pix_get_depth(pix) == 1
        true
    else
        false
    end
end

# =================================================================================================
"""
    pix_write_tiff_i(
        pix::Pix,
        compression::IFF
    )::Tuple{Ptr{UInt8}, Csize_t}

Write an image to memory in the PNG image format. If there is an error `(C_NULL, 0)` is returned.

__Arguments:__

| T | Name        | Default  | Description
|:--| :---------- | :------- | :----------
| R | pix         |          | The image to write to memory.
| R | compression |          | The compression to use on the image.

__Restrictions:__

  * `compression` - Must be a one of the following compression formats:

      * `IFF_TIFF` - Supports all images.
      * `IFF_TIFF_RLE` - Requires a b&w 1bpp image.
      * `IFF_TIFF_PACKBITS` - Requires a b&w 1bpp image.
      * `IFF_TIFF_G3` - Requires a b&w 1bpp image.
      * `IFF_TIFF_G4` - Requires a b&w 1bpp image.
      * `IFF_TIFF_LZW` - Supports all images.
      * `IFF_TIFF_ZIP` - Supports all images.
      * `IFF_TIFF_JPEG` - Supports all images.

__Details:__

**This method is meant for internal use only.**  It returns the memory pointer and number of bytes
allocated by the Leptonica library to write out the image.  The caller is responsible for freeing
this memory with a call to `lept_free()`.
"""
function pix_write_tiff_i(
            pix::Pix,
            compression::IFF
        )::Tuple{Ptr{UInt8}, Csize_t}
    local temp = Ref(Ptr{UInt8}(C_NULL))
    local size = Ref(Csize_t(0))

    # ---------------------------------------------------------------------------------------------
    # Verify the input values.
    if is_valid(pix) == false
        @error "Pix has been freed."
        return (C_NULL, 0)
    end

    if compression ∉ TIFF_FORMATS
        @error "Invalid compression format."
        return (C_NULL, 0)
    end

    if valid_tiff_compression(pix, compression) == false
        @error "Invalid compression, bit depth must be 1 for compression $compression."
        return (C_NULL, 0)
    end

    # ---------------------------------------------------------------------------------------------
    # Perform the write.
    local retval = ccall(
        (:pixWriteMemTiff, LEPTONICA),
        Cint,
        (Ptr{Ptr{UInt8}}, Ptr{Csize_t}, Ptr{Cvoid}, Cint),
        temp,
        size,
        pix,
        compression
    )

    if retval == 1
        return (C_NULL, 0)
    end

    return (temp[], size[])
end

# =================================================================================================
"""
    pix_write_tiff(
        filename::AbstractString,
        pix::Pix;
        compression::IFF = IFF_TIFF,
        append::Bool = false
    )::Bool

Write an image to a file in the TIFF image format. Returns `false` if there was an error.

__Arguments:__

| T | Name        | Default    | Description
|:--| :---------- | :--------- | :----------
| R | filename    |            | The name of the file to write to or create.
| R | pix         |            | The image to write to disk.
| O | compression | `IFF_TIFF` | The compression to use on the image.
| O | append      | `false`    | Should we overwrite the image or append it to the file.

__Restrictions:__

  * `compression` - Must be a one of the following compression formats:

      * `IFF_TIFF` - Supports all images.
      * `IFF_TIFF_RLE` - Requires a b&w 1bpp image.
      * `IFF_TIFF_PACKBITS` - Requires a b&w 1bpp image.
      * `IFF_TIFF_G3` - Requires a b&w 1bpp image.
      * `IFF_TIFF_G4` - Requires a b&w 1bpp image.
      * `IFF_TIFF_LZW` - Supports all images.
      * `IFF_TIFF_ZIP` - Supports all images.
      * `IFF_TIFF_JPEG` - Supports all images.

__Details:__

The default compression is no compression.

TIFF files can contain multiple image.  If a file does not exist it will be created.  If the file
exists and append is `false` then the file will be overwritten.  If append is `true` then the image
will be added to the file.
"""
function pix_write_tiff(
            filename::AbstractString,
            pix::Pix;
            compression::IFF = IFF_TIFF,
            append::Bool = false
        )::Bool

    if is_valid(pix) == false
        @error "Pix has been freed."
        return false
    end

    if compression ∉ TIFF_FORMATS
        @error "Invalid compression format."
        return false
    end

    if valid_tiff_compression(pix, compression) == false
        @error "Invalid compression, bit depth must be 1 for compression $compression."
        return false
    end

    local retval = ccall(
        (:pixWriteTiff, LEPTONICA),
        Cint,
        (Cstring, Ptr{Cvoid}, Cint, Cstring),
        filename,
        pix,
        compression,
        append ? "a" : "w"
    )

    return retval == 0
end

# =================================================================================================
"""
    pix_write_tiff(
        stream::IO,
        pix::Pix;
        compression::IFF = IFF_TIFF
    )::Bool

Write an image to an IO stream in the TIFF image format. Returns `false` if there was an error.

__Arguments:__

| T | Name        | Default    | Description
|:--| :---------- | :--------- | :----------
| R | stream      |            | The stream to write the image to.
| R | pix         |            | The image to write to disk.
| O | compression | `IFF_TIFF` | The compression to use on the image.

__Restrictions:__

  * `compression` - Must be a one of the following compression formats:

      * `IFF_TIFF` - Supports all images.
      * `IFF_TIFF_RLE` - Requires a b&w 1bpp image.
      * `IFF_TIFF_PACKBITS` - Requires a b&w 1bpp image.
      * `IFF_TIFF_G3` - Requires a b&w 1bpp image.
      * `IFF_TIFF_G4` - Requires a b&w 1bpp image.
      * `IFF_TIFF_LZW` - Supports all images.
      * `IFF_TIFF_ZIP` - Supports all images.
      * `IFF_TIFF_JPEG` - Supports all images.

__Details:__

The default compression is no compression.
"""
function pix_write_tiff(
            stream::IO,
            pix::Pix;
            compression::IFF = IFF_TIFF
        )::Bool
    local result = false
    local data, size = pix_write_tiff_i(pix, compression)

    if data != C_NULL
        try
            if size > 0
                unsafe_write(stream, data, size)
                result = true
            end
        finally
            lept_free(data)
        end
    end

    return result
end

# =================================================================================================
"""
    pix_write_tiff(
        pix::Pix;
        compression::IFF = IFF_TIFF
    )::Union{Vector{UInt8}, Nothing}

Write an image to an byte array in the TIFF image format. Returns `nothing` if there was an error.

__Arguments:__

| T | Name        | Default    | Description
|:--| :---------- | :--------- | :----------
| R | pix         |            | The image to write to a byte array.
| O | compression | `IFF_TIFF` | The compression to use on the image.

__Restrictions:__

  * `compression` - Must be a one of the following compression formats:

      * `IFF_TIFF` - Supports all images.
      * `IFF_TIFF_RLE` - Requires a b&w 1bpp image.
      * `IFF_TIFF_PACKBITS` - Requires a b&w 1bpp image.
      * `IFF_TIFF_G3` - Requires a b&w 1bpp image.
      * `IFF_TIFF_G4` - Requires a b&w 1bpp image.
      * `IFF_TIFF_LZW` - Supports all images.
      * `IFF_TIFF_ZIP` - Supports all images.
      * `IFF_TIFF_JPEG` - Supports all images.

__Details:__

The default compression is no compression.
"""
function pix_write_tiff(
            pix::Pix;
            compression::IFF = IFF_TIFF
        )::Union{Vector{UInt8}, Nothing}
    local output     = nothing
    local data, size = pix_write_tiff_i(pix, compression)

    if data != C_NULL
        try
            if size > 0
                output = Vector{UInt8}(undef, size)
                unsafe_copyto!(pointer(output), data, size)
            end
        finally
            lept_free(data)
        end
    end

    return output
end


text:
function read_integers()
    [parse(Int, x) for x in split(readline())]
end

function read_instance()
    n, m = read_integers()
    conditions = []
    for i=1:m
        push!(conditions, read_integers())
    end
    k = read_integers()[1]
    choices = []
    for i=1:k
        push!(choices, read_integers())
    end
    n, m, conditions, k, choices
end

function put_ball(choices, x, n, k)
    balls = zeros(Int, n)
    for i=0:k-1
        if x >> i & 1 == 1
            balls[choices[i+1][1]] += 1
        else
            balls[choices[i+1][2]] += 1
        end
    end
    balls
end

function num_satisfied(conditions, balls)
    ret = 0
    for cond in conditions
        if balls[cond[1]] ≥ 1 && balls[cond[2]] ≥ 1
            ret += 1
        end
    end
    ret
end

function solve()
    n, m, conditions, k, choices = read_instance()
    ret = 0
    for x=0:2^k
        balls = put_ball(choices, x, n, k)
        ret = max(ret, num_satisfied(conditions, balls))
    end
    ret
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


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

using Markdown
using InteractiveUtils

# ╔═╡ 571c959e-1aee-11eb-0e1e-a1b4d8f7762d
using Distributions, LaTeXStrings, Plots

# ╔═╡ f7c7716c-3593-11eb-00c3-677a7dbc5fec
md"""
# Plots comparing bounds
"""

# ╔═╡ 6e5aab50-b74b-11eb-33c3-5d53afe92e95
md"""
## Problem 1: any number of stages, any decoder, etc.
"""

# ╔═╡ 7fdcfedc-b74b-11eb-399f-95f6fdf61e67
md"""
### Lower bound (binary entropy)
"""

# ╔═╡ 995ee74e-b74b-11eb-20ef-0fac0d588eac
prob1_lb(p) = entropy(Bernoulli(p),2)

# ╔═╡ fb810c18-b74b-11eb-39a6-5b769158670f
md"""
### Upper bound for fully adaptive

Theorem 3 from:
Aldridge, Matthew. "Rates of adaptive group testing in the linear regime." 2019 IEEE International Symposium on Information Theory (ISIT). IEEE, 2019.
"""

# ╔═╡ bd23cdb0-b74c-11eb-1d4e-2bd846b6d795
function adap_ub(p)
	q = 1-p

	m = ceil(-log(2-p)/log(1-p))
	a = floor(log(2,m))
	b = m - 2^a

	F = q^m + (1-q^(m-2*b))*(a+1) + (q^(m-2*b)-q^m)*(a+2)
	G = m*q^m + 1/p*(1 + m*q^(m+1) - (m+1)*q^m)
	return F/G
end

# ╔═╡ aac6b044-b74a-11eb-2199-1b6613ccaa1b
md"""
## Problem 2: two-stage conservative testing
"""

# ╔═╡ 1ac80b88-3594-11eb-2a77-7321f527f86a
md"""
### Lower bound

Theorem 5 from:
Aldridge, Matthew. "Conservative two-stage group testing." arXiv preprint arXiv:2005.06617 (2020).
"""

# ╔═╡ 0919ebde-b74b-11eb-3947-878a1f65f615
function prob2_lb(p;wrange=2:100000)
	if p >= (3-sqrt(5))/2
		return 1.0
	else
		fp = maximum(w->-w*log(1-(1-p)^(w-1)),wrange)
		gp = maximum(w->-w*log(1-(1-p)^w    ),wrange)
		bound2 = 1/gp*(log(gp)+1)
		bound3 = p + 1/fp*(log((1-p)*fp)+1)
		return max(bound2,bound3)
	end
end

# ╔═╡ 95bfab7c-b74a-11eb-04d7-ab9d03e13c5b
md"""
### Upper bound for HYPER
"""

# ╔═╡ 9e393e90-35f0-11eb-25bc-af3fc684f85f
Base.binomial(m::Float64,q::Int) = q == 0 ? 1 : prod(m-i for i in 0:q-1)/factorial(q)

# ╔═╡ 1db68090-3594-11eb-3c00-1f948011570b
function hyper_ub(p,n,m,q)
	α = 0
	β = 1
	k = n*q/m
	r = 1-p
	u = binomial(m-2,q-2)*ceil(n/binomial(m,q))
	p1 = 1 - β + (β-α)*r^k
	p2 = p1^2 + (β-α)^2*r^(2*k-u)*(1-r^u)
	return (m + n*(1 - q*p1 + binomial(q,2)*p2))/n
end

# ╔═╡ 215f9960-3595-11eb-313d-ef186d77265b
md"""
Numerical optimization of bound:
"""

# ╔═╡ 5c9a1260-3595-11eb-0e60-3fe0c174ca5a
hyper_ub(p,n,q;maxm=n) = minimum(hyper_ub(p,n,m,q) for m in q:maxm)

# ╔═╡ 2e3ff118-b74b-11eb-2204-034129ec8abe
hyper_ub(p,n;qrange=1:3) = minimum(hyper_ub(p,n,q) for q in qrange)

# ╔═╡ 2383c37a-35a8-11eb-3f37-ed1672ed16dc
md"""
## Plot
"""

# ╔═╡ 44682866-b1ee-11eb-15dc-69a15ac82723
let n = 6144
	pal = palette([RGB(0,0,0); PlotThemes.wong_palette])
	prange = 0.001:0.001:0.2
	pticks = 0.01:0.01:0.2
	plot(size=(800,400),dpi=200,legend=:bottomright,
		xlabel=L"\textrm{Prevalence } (\%)",
		ylabel=L"\textrm{tests/individual}: \mathrm{E}\;T/n",
		xticks=(pticks,string.(convert.(Int,pticks*100))),xlims=(0,maximum(pticks)),
		yticks=0:0.2:1,ylims=(0,1),
		palette=pal
	)

	# Problem 2
	plot!(prange,p->hyper_ub(p,n);label="Upper bound for HYPER method",linewidth=3)
	plot!(prange,prob2_lb;label="Lower bound for problem 2: two-stage testing with conservative decoder",linewidth=3,linestyle=:dash)

	# Problem 1
	plot!(prange,adap_ub;label="Upper bound for fully adaptive method",linewidth=3)
	plot!(prange,prob1_lb;label="Lower bound for problem 1: any number of stages, any decoder, etc.",linewidth=3,linestyle=:dash)

	savefig("fig.png"); plot!()
end

# ╔═╡ Cell order:
# ╟─f7c7716c-3593-11eb-00c3-677a7dbc5fec
# ╠═571c959e-1aee-11eb-0e1e-a1b4d8f7762d
# ╟─6e5aab50-b74b-11eb-33c3-5d53afe92e95
# ╟─7fdcfedc-b74b-11eb-399f-95f6fdf61e67
# ╠═995ee74e-b74b-11eb-20ef-0fac0d588eac
# ╟─fb810c18-b74b-11eb-39a6-5b769158670f
# ╠═bd23cdb0-b74c-11eb-1d4e-2bd846b6d795
# ╟─aac6b044-b74a-11eb-2199-1b6613ccaa1b
# ╟─1ac80b88-3594-11eb-2a77-7321f527f86a
# ╠═0919ebde-b74b-11eb-3947-878a1f65f615
# ╟─95bfab7c-b74a-11eb-04d7-ab9d03e13c5b
# ╠═9e393e90-35f0-11eb-25bc-af3fc684f85f
# ╠═1db68090-3594-11eb-3c00-1f948011570b
# ╟─215f9960-3595-11eb-313d-ef186d77265b
# ╠═5c9a1260-3595-11eb-0e60-3fe0c174ca5a
# ╠═2e3ff118-b74b-11eb-2204-034129ec8abe
# ╟─2383c37a-35a8-11eb-3f37-ed1672ed16dc
# ╠═44682866-b1ee-11eb-15dc-69a15ac82723


text:
function minimize!(objective::ObjWrapper, x0,
                   scheme::QuasiNewton, options=OptOptions(),
                   cache=preallocate_qn_caches_inplace(x0))

    minimize!(objective, (x0, nothing), (scheme, Backtracking()), options, cache)
end
function minimize!(objective::ObjWrapper, x0,
                   approach::Tuple, options=OptOptions(),
                   cache=preallocate_qn_caches_inplace(x0))

    minimize!(objective, (x0, nothing), approach, options, cache)
end
function minimize!(objective::ObjWrapper, state0::Tuple,
                   scheme::QuasiNewton, options=OptOptions(),
                   cache=preallocate_qn_caches_inplace(first(state0)))

    minimize!(objective, state0, (scheme, Backtracking()), options, cache)
end
function minimize!(objective::ObjWrapper, state0::Tuple, approach::Tuple{<:Any, <:LineSearch},
                   options::OptOptions=OptOptions(), # options
                   cache = preallocate_qn_caches_inplace(x0), # preallocate arrays for QN
                   )

    #==============
         Setup
    ==============#
    x0, B0 = state0
    T = eltype(x0)
    x, fx, ∇fx, z, fz, ∇fz, B = prepare_variables(objective, approach, x0, copy(x0), B0)

    #========================
         First iteration
    ========================#
    x, fx, ∇fx, z, fz, ∇fz, B, is_converged = iterate!(cache, x, fx, ∇fx, z, fz, ∇fz, B, approach, objective, options)

    iter = 0
    while iter <= options.maxiter && !is_converged
        iter += 1
        # take a step and update approximation
        x, fx, ∇fx, z, fz, ∇fz, B, is_converged = iterate!(cache, x, fx, ∇fx, z, fz, ∇fz, B, approach, objective, options, false)
        if norm(x.-z, Inf) ≤ eps(T)
            break
        end
    end
    return z, fz, ∇fz, iter
end

function iterate!(cache, x, fx::Tf, ∇fx, z, fz, ∇fz, B, approach::Tuple{<:Any, <:LineSearch}, objective, options, is_first=nothing) where Tf
    # split up the approach into the hessian approximation scheme and line search
    scheme, linesearch = approach

    y, d, s = cache.y, cache.d, cache.s

    # Move nexts into currs
    fx = fz
    copyto!(x, z)
    copyto!(∇fx, ∇fz)

    # Update current gradient and calculate the search direction
    d = find_direction!(d, B, ∇fx, scheme) # solve Bd = -∇fx
    φ = LineObjective!(objective, ∇fz, z, x, d, fx, dot(∇fx, d))

    # Perform line search along d
    α, f_α, ls_success = find_steplength(linesearch, φ, Tf(1.0))

    # Calculate final step vector and update the state
    @. s = α * d
    @. z = x + s

    # Update approximation
    fz, ∇fz, B = update_obj!(objective, d, s, y, ∇fx, z, ∇fz, B, scheme, is_first)

    # Check for gradient convergence
    is_converged = converged(z, ∇fz, options.g_tol)

    return x, fx, ∇fx, z, fz, ∇fz, B, is_converged
end


text:
include("../src/JOH.jl")
include("../GAP/GAP.jl")

GAP.load_folder()

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
module SparseBenchmarks

include(joinpath(dirname(@__FILE__), "..", "utils", "RandUtils.jl"))

using .RandUtils
using BenchmarkTools
using SparseArrays
using LinearAlgebra

const SUITE = BenchmarkGroup(["array"])

#########
# index #
#########

# Note that some of the "logical" tests are commented
# out because they require resolution of JuliaLang/julia#14717.

g = addgroup!(SUITE, "index")

# vector #
#--------#

sizes = (10^3, 10^4, 10^5)

spvecs = map(s -> samesprand(s, inv(sqrt(s))), sizes)
splogvecs = map(s -> samesprandbool(s, 1e-5), sizes)

for (s, v, l) in zip(sizes, spvecs, splogvecs)
    g["spvec", "array",   s] = @benchmarkable getindex($v, $(samerand(1:s, s)))
    g["spvec", "integer", s] = @benchmarkable getindex($v, $(samerand(1:s)))
    g["spvec", "range",   s] = @benchmarkable getindex($v, $(1:s))
    g["spvec", "logical", s] = @benchmarkable getindex($v, $(samerand(Bool, s)))
    # g["spvec", "splogical", s, nnz(v), nnz(l)] = @benchmarkable getindex($v, $l)
end

# matrix #
#--------#

sizes = (10, 10^2, 10^3)
inds = map(s -> samerand(1:s), sizes)
matrices = map(s -> samesprand(s, s, inv(sqrt(s))), sizes)
vectors = map(s -> samerand(1:s, s), sizes)
logvecs = map(s -> samerand(Bool, s), sizes)
splogmats = map(s -> samesprandbool(s, s, 1e-5), sizes)
splogvecs = map(s -> samesprandbool(s, 1, 1e-5), sizes)

for (s, m, v, l, sl, c) in zip(sizes, matrices, vectors, logvecs, splogvecs, inds)
    g["spmat", "col", "array", s] = @benchmarkable getindex($m, $v, $c)
    g["spmat", "col", "range", s] = @benchmarkable getindex($m, $(1:s), $c)
    g["spmat", "col", "OneTo", s] = @benchmarkable getindex($m, $(Base.OneTo(s)), $c)
    g["spmat", "col", "logical", s] = @benchmarkable getindex($m, $l, $c)
    # g["spmat", "col", "splogical", s] = @benchmarkable getindex($m, $sl, $c)
end

for (s, m, v, l, sl, r) in zip(sizes, matrices, vectors, logvecs, splogvecs, inds)
    g["spmat", "row", "array", s] = @benchmarkable getindex($m, $r, $v)
    g["spmat", "row", "range", s] = @benchmarkable getindex($m, $r, $(1:s))
    g["spmat", "row", "OneTo", s] = @benchmarkable getindex($m, $r, $(Base.OneTo(s)))
    g["spmat", "row", "logical", s] = @benchmarkable getindex($m, $r, $l)
    # g["spmat", "row", "splogical", s] = @benchmarkable getindex($m, $r, $sl)
end

for (s, m, v, l, sl, i) in zip(sizes, matrices, vectors, logvecs, splogmats, inds)
    g["spmat", "array", s] = @benchmarkable getindex($m, $v, $v)
    g["spmat", "integer", s] = @benchmarkable getindex($m, $i, $i)
    g["spmat", "range", s] = @benchmarkable getindex($m, $(1:s), $(1:s))
    g["spmat", "OneTo", s] = @benchmarkable getindex($m, $(Base.OneTo(s)), $(Base.OneTo(s)))
    g["spmat", "logical", s] = @benchmarkable getindex($m, $l, $l)
    g["spmat", "splogical", s] = @benchmarkable getindex($m, $sl)
end

for b in values(g)
    b.params.time_tolerance = 0.3
end

######################
# transpose (#14631) #
######################

small_sqr = samesprand(600, 600, 0.01)
small_rct = samesprand(600, 400, 0.01)
large_sqr = samesprand(20000, 20000, 0.01)
large_rct = samesprand(20000, 10000, 0.01)

g = addgroup!(SUITE, "transpose", ["adjoint"])

for m in (small_sqr, small_rct, large_sqr, large_rct)
    cm = m + m*im
    s = size(m)
    g["transpose", s] = @benchmarkable transpose($m)
    g["transpose!", s] = @benchmarkable transpose!($(copy(transpose(m))), $m)
    g["adjoint", s] = @benchmarkable adjoint($cm)
    g["adjoint!", s] = @benchmarkable adjoint!($(copy(transpose(cm))), $cm)
end

for b in values(g)
    b.params.time_tolerance = 0.3
end

##############
# arithmetic #
##############

g = addgroup!(SUITE, "arithmetic")

# unary minus, julialang repo issue #19503 / fix #19530
g["unary minus", size(small_sqr)] = @benchmarkable -$small_sqr
g["unary minus", size(large_sqr)] = @benchmarkable -$large_sqr

for b in values(g)
    b.params.time_tolerance = 0.3
end

################
# constructors #
################
g = addgroup!(SUITE, "constructors")

const UPLO = :U
for s in sizes
    nz = floor(Int, 1e-4*s*s)
    I = samerand(1:s, nz)
    J = samerand(1:s, nz)
    V = randvec(nz)
    g["IV", s] = @benchmarkable sparsevec($I, $V)
    g["IJV", s] = @benchmarkable sparse($I, $J, $V)
    g["Diagonal", s] = @benchmarkable sparse($(Diagonal(randvec(s))))
    g["Bidiagonal", s] = @benchmarkable sparse($(Bidiagonal(randvec(s), randvec(s-1), UPLO)))
    g["Tridiagonal", s] = @benchmarkable sparse($(Tridiagonal(randvec(s-1), randvec(s), randvec(s-1))))
    g["SymTridiagonal", s] = @benchmarkable sparse($(SymTridiagonal(randvec(s), randvec(s-1))))
end

#########################
# matrix multiplication #
#########################

g = addgroup!(SUITE, "matmul")

# mixed sparse-dense matmul #
#---------------------------#

using LinearAlgebra: *, mul!

function allocmats_ds(om, ok, on, s, nnzc, T)
    m, k, n = map(x -> Int(s*x), (om, ok, on))
    densemat, sparsemat = samerand(T, m, k), samesprand(T, k, n, nnzc/k)
    tdensemat = transpose!(similar(densemat, reverse(size(densemat))), densemat)
    tsparsemat = transpose!(similar(sparsemat, reverse(size(sparsemat))), sparsemat)
    destmat = similar(densemat, m, n)
    return m, k, n, destmat,
        densemat, sparsemat,
        tdensemat, tsparsemat
end
function allocmats_sd(om, ok, on, s, nnzc, T)
    m, k, n = map(x -> Int(s*x), (om, ok, on))
    densemat, sparsemat = samerand(T, k, m), samesprand(T, n, k, nnzc/n)
    tdensemat = transpose!(similar(densemat, reverse(size(densemat))), densemat)
    tsparsemat = transpose!(similar(sparsemat, reverse(size(sparsemat))), sparsemat)
    destmat = similar(densemat, n, m)
    return m, k, n, destmat,
        densemat, sparsemat,
        tdensemat, tsparsemat
end

for (om, ok, on) in (# order of matmul dimensions m, k, and n
        (10^2, 10^2, 10^2),  # dense square * sparse square -> dense square
        (10^1, 10^1, 10^3),  # dense square * sparse short -> dense short
        (10^2, 10^2, 10^1),  # dense square * sparse tall -> dense tall
        (10^1, 10^3, 10^3),  # dense short * sparse square -> dense short
        (10^1, 10^2, 10^3),  # dense short * sparse short -> dense short
        (10^1, 10^3, 10^2),  # dense short * sparse tall -> dense short
        (10^3, 10^1, 10^1),  # dense tall * sparse square -> dense tall
        (10^2, 10^1, 10^2),  # dense tall * sparse short -> dense square
        ) # the preceding descriptions apply to dense-sparse matmul without
          # any transpositions. other cases are described below
    #
    # the transpose and adjoint variants share kernel code
    # the in-place and out-of-place variants share kernel code
    # so exercise the different variants in different ways
    #
    # for A[t|c]_mul_B[t|c][!]([dense,], dense, sparse) kernels,
    # the dense matrix is m-by-k, or k-by-m for A(c|t) operations
    # the sparse matrix is k-by-n, or n-by-k for B(c|t) operations
    # and the (dense) destination matrix is m-by-n in any case
    # the sparse matrix has approximately 10 entries per column
    #
    # # out-of-place dense-sparse ops, transpose variants, i.e. A[t]_mul_B[t](dense, sparse)
    m, k, n, destmat, densemat, sparsemat, tdensemat, tsparsemat = allocmats_ds(om, ok, on, 1/2, 4, Float64)
    g["A_mul_B",   "dense $(m)x$(k), sparse $(k)x$(n) -> dense $(m)x$(n)"] = @benchmarkable *($densemat, $sparsemat)
    g["A_mul_Bt",  "dense $(m)x$(k), sparse $(n)x$(k) -> dense $(m)x$(n)"] = @benchmarkable *($densemat, $(Transpose(tsparsemat)))
    g["At_mul_B",  "dense $(k)x$(m), sparse $(k)x$(n) -> dense $(m)x$(n)"] = @benchmarkable *($(Transpose(tdensemat)), $sparsemat)
    g["At_mul_Bt", "dense $(k)x$(m), sparse $(n)x$(k) -> dense $(m)x$(n)"] = @benchmarkable *($(Transpose(tdensemat)), $(Transpose(tsparsemat)))
    # in-place dense-sparse -> dense ops, transpose variants, i.e. A[t]_mul[t]!(dense, dense, sparse)
    m, k, n, destmat, densemat, sparsemat, tdensemat, tsparsemat = allocmats_ds(om, ok, on, 4, 12, Float64)
    g["A_mul_B!",   "dense $(m)x$(k), sparse $(k)x$(n) -> dense $(m)x$(n)"] = @benchmarkable mul!($destmat, $densemat, $sparsemat)
    g["A_mul_Bt!",  "dense $(m)x$(k), sparse $(n)x$(k) -> dense $(m)x$(n)"] = @benchmarkable mul!($destmat, $densemat, $(Transpose(tsparsemat)))
    g["At_mul_B!",  "dense $(k)x$(m), sparse $(k)x$(n) -> dense $(m)x$(n)"] = @benchmarkable mul!($destmat, $(Transpose(tdensemat)), $sparsemat)
    g["At_mul_Bt!", "dense $(k)x$(m), sparse $(n)x$(k) -> dense $(m)x$(n)"] = @benchmarkable mul!($destmat, $(Transpose(tdensemat)), $(Transpose(tsparsemat)))
    # out-of-place dense-sparse ops, adjoint variants, i.e. A[c]_mul_B[c](dense, sparse)
    m, k, n, destmat, densemat, sparsemat, tdensemat, tsparsemat = allocmats_ds(om, ok, on, 1/2, 4, Complex{Float64})
    g["A_mul_Bc",  "dense $(m)x$(k), sparse $(n)x$(k) -> dense $(m)x$(n)"] = @benchmarkable *($densemat, $(Adjoint(tsparsemat)))
    g["Ac_mul_B",  "dense $(k)x$(m), sparse $(k)x$(n) -> dense $(m)x$(n)"] = @benchmarkable *($(Adjoint(tdensemat)), $sparsemat)
    g["Ac_mul_Bc", "dense $(k)x$(m), sparse $(n)x$(k) -> dense $(m)x$(n)"] = @benchmarkable *($(Adjoint(tdensemat)), $(Adjoint(tsparsemat)))
    # in-place dense-sparse -> dense ops, adjoint variants, i.e. A[c]_mul[c]!(dense, dense, sparse)
    m, k, n, destmat, densemat, sparsemat, tdensemat, tsparsemat = allocmats_ds(om, ok, on, 2, 8, Complex{Float64})
    g["A_mul_Bc!",  "dense $(m)x$(k), sparse $(n)x$(k) -> dense $(m)x$(n)"] = @benchmarkable mul!($destmat, $densemat, $(Adjoint(tsparsemat)))
    g["Ac_mul_B!",  "dense $(k)x$(m), sparse $(k)x$(n) -> dense $(m)x$(n)"] = @benchmarkable mul!($destmat, $(Adjoint(tdensemat)), $sparsemat)
    g["Ac_mul_Bc!", "dense $(k)x$(m), sparse $(n)x$(k) -> dense $(m)x$(n)"] = @benchmarkable mul!($destmat, $(Adjoint(tdensemat)), $(Adjoint(tsparsemat)))
    #
    # for A[t|c]_mul_B[t|c][!]([dense,], sparse, dense) kernels,
    # the sparse matrix is n-by-k, or k-by-n for B(c|t) operations
    # the dense matrix is k-by-m, or m-by-k for A(c|t) operations
    # and the (dense) destination matrix is n-by-m in any case
    # the sparse matrix has approximately 10 entries per column
    #
    # out-of-place sparse-dense ops, transpose variants, i.e. A[t]_mul_B[t](sparse, dense)
    m, k, n, destmat, densemat, sparsemat, tdensemat, tsparsemat = allocmats_sd(om, ok, on, 1/2, 4, Complex{Float64})
    g["A_mul_B",   "sparse $(n)x$(k), dense $(k)x$(m) -> dense $(n)x$(m)"] = @benchmarkable *($sparsemat, $densemat)
    g["A_mul_Bt",  "sparse $(n)x$(k), dense $(m)x$(k) -> dense $(n)x$(m)"] = @benchmarkable *($sparsemat, $(Transpose(tdensemat)))
    g["At_mul_B",  "sparse $(k)x$(n), dense $(k)x$(m) -> dense $(n)x$(m)"] = @benchmarkable *($(Transpose(tsparsemat)), $densemat)
    g["At_mul_Bt", "sparse $(k)x$(n), dense $(m)x$(k) -> dense $(n)x$(m)"] = @benchmarkable *($(Transpose(tsparsemat)), $(Transpose(tdensemat)))
    # in-place sparse-dense -> dense ops, transpose variants, i.e. A[t|c]_mul_B[t|c]!(dense, sparse, dense)
    m, k, n, destmat, densemat, sparsemat, tdensemat, tsparsemat = allocmats_sd(om, ok, on, 4, 12, Complex{Float64})
    g["A_mul_B!",   "sparse $(n)x$(k), dense $(k)x$(m) -> dense $(n)x$(m)"] = @benchmarkable mul!($destmat, $sparsemat, $densemat)
    g["A_mul_Bt!",  "sparse $(n)x$(k), dense $(m)x$(k) -> dense $(n)x$(m)"] = @benchmarkable mul!($destmat, $sparsemat, $(Transpose(tdensemat)))
    g["At_mul_B!",  "sparse $(k)x$(n), dense $(k)x$(m) -> dense $(n)x$(m)"] = @benchmarkable mul!($destmat, $(Transpose(tsparsemat)), $densemat)
    g["At_mul_Bt!", "sparse $(k)x$(n), dense $(m)x$(k) -> dense $(n)x$(m)"] = @benchmarkable mul!($destmat, $(Transpose(tsparsemat)), $(Transpose(tdensemat)))
    # out-of-place sparse-dense ops, adjoint variants, i.e. A[c]_mul_B[c](sparse, dense)
    m, k, n, destmat, densemat, sparsemat, tdensemat, tsparsemat = allocmats_sd(om, ok, on, 1/2, 4, Complex{Float64})
    g["A_mul_Bc",  "sparse $(n)x$(k), dense $(m)x$(k) -> dense $(n)x$(m)"] = @benchmarkable *($sparsemat, $(Adjoint(tdensemat)))
    g["Ac_mul_B",  "sparse $(k)x$(n), dense $(k)x$(m) -> dense $(n)x$(m)"] = @benchmarkable *($(Adjoint(tsparsemat)), $densemat)
    g["Ac_mul_Bc", "sparse $(k)x$(n), dense $(m)x$(k) -> dense $(n)x$(m)"] = @benchmarkable *($(Adjoint(tsparsemat)), $(Adjoint(tdensemat)))
    # in-place sparse-dense -> dense ops, adjoint variants, i.e. A[t|c]_mul_B[t|c]!(dense, sparse, dense)
    m, k, n, destmat, densemat, sparsemat, tdensemat, tsparsemat = allocmats_sd(om, ok, on, 2, 8, Complex{Float64})
    g["A_mul_Bc!",  "sparse $(n)x$(k), dense $(m)x$(k) -> dense $(n)x$(m)"] = @benchmarkable mul!($destmat, $sparsemat, $(Adjoint(tdensemat)))
    g["Ac_mul_B!",  "sparse $(k)x$(n), dense $(k)x$(m) -> dense $(n)x$(m)"] = @benchmarkable mul!($destmat, $(Adjoint(tsparsemat)), $densemat)
    g["Ac_mul_Bc!", "sparse $(k)x$(n), dense $(m)x$(k) -> dense $(n)x$(m)"] = @benchmarkable mul!($destmat, $(Adjoint(tsparsemat)), $(Adjoint(tdensemat)))
end

for b in values(g)
    b.params.time_tolerance = 0.3
end

#################
# sparse matvec #
#################
g = addgroup!(SUITE, "sparse matvec")
B = randn(100000, 100)
A = sprand(100000, 100000, 0.00001)
g["non-adjoint"] = @benchmarkable A * B
g["adjoint"] = @benchmarkable A' * B

#################
# sparse solves #
#################
g = addgroup!(SUITE, "sparse solves")
# Problem similar to issue #30288
let m = 10000, n = 9000
    A = spdiagm(0 => fill(2.0, m),
               -1 => fill(1.0, m - 1),
                1 => fill(1.0, m - 1),
              360 => fill(1.0, m - 360))[:, 1:n]
    AtA = A'A
    b   = ones(m)
    B   = ones(m, 2)
    Atb = A'b
    AtB = A'B

    g["least squares (default), vector rhs"]  = @benchmarkable         $A\$b
    g["least squares (default), matrix rhs"]  = @benchmarkable         $A\$B
    g["least squares (qr), vector rhs"]   = @benchmarkable     qr($A)\$b
    g["least squares (qr), matrix rhs"]   = @benchmarkable     qr($A)\$B
    g["square system (default), vector rhs"]  = @benchmarkable       $AtA\$Atb
    g["square system (default), matrix rhs"]  = @benchmarkable       $AtA\$AtB
    g["square system (ldlt), vector rhs"] = @benchmarkable ldlt($AtA)\$Atb
    g["square system (ldlt), matrix rhs"] = @benchmarkable ldlt($AtA)\$AtB
    g["square system (lu), vector rhs"]   = @benchmarkable   lu($AtA)\$Atb
    g["square system (lu), matrix rhs"]   = @benchmarkable   lu($AtA)\$AtB
end

end # module


text:
module Breakpoints

using ..JuliaInterpreter
using JuliaInterpreter: JuliaFrameCode, JuliaStackFrame, BreakpointState,
                        truecondition, falsecondition, prepare_framecode, get_framecode,
                        sparam_syms, linenumber, statementnumber
using Base.Meta: isexpr
using InteractiveUtils

export @breakpoint, breakpoint, enable, disable, remove

# A type that is unique to this package for which there are no valid operations
struct Unassigned end

"""
    BreakpointRef(framecode, stmtidx)
    BreakpointRef(framecode, stmtidx, err)

A reference to a breakpoint at a particular statement index `stmtidx` in `framecode`.
If the break was due to an error, supply that as well.
"""
struct BreakpointRef
    framecode::JuliaFrameCode
    stmtidx::Int
    err
end
BreakpointRef(framecode, stmtidx) = BreakpointRef(framecode, stmtidx, nothing)

function Base.show(io::IO, bp::BreakpointRef)
    if checkbounds(Bool, bp.framecode.breakpoints, bp.stmtidx)
        lineno = linenumber(bp.framecode, bp.stmtidx)
        print(io, "breakpoint(", bp.framecode.scope, ", ", lineno)
    else
        print(io, "breakpoint(", bp.framecode.scope, ", %", bp.stmtidx)
    end
    if bp.err !== nothing
        print(io, ", ", bp.err)
    end
    print(io, ')')
end

const _breakpoints = BreakpointRef[]
breakpoints() = copy(_breakpoints)

Base.getindex(bp::BreakpointRef) = bp.framecode.breakpoints[bp.stmtidx]
function Base.setindex!(bp::BreakpointRef, isactive::Bool)
    bp.framecode.breakpoints[bp.stmtidx] = BreakpointState(isactive, bp[].condition)
end
function toggle!(bp::BreakpointRef)
    state = bp[]
    bp.framecode.breakpoints[bp.stmtidx] = BreakpointState(!state.isactive, state.condition)
end

function add_breakpoint(framecode, stmtidx)
    bp = BreakpointRef(framecode, stmtidx)
    # Since there can be only one BreakpointState for a given framecode/stmtidx,
    # check whether _breakpoints is already storing a reference to that location
    idx = findfirst(isequal(bp), _breakpoints)
    if idx === nothing
        push!(_breakpoints, bp)
    end
    return bp
end

function shouldbreak(frame, pc=frame.pc[])
    idx = convert(Int, pc)
    isassigned(frame.code.breakpoints, idx) || return false
    bp = frame.code.breakpoints[idx]
    bp.isactive || return false
    return Base.invokelatest(bp.condition, frame)::Bool
end

function prepare_slotfunction(framecode::JuliaFrameCode, body::Union{Symbol,Expr})
    ismeth = framecode.scope isa Method
    uslotnames = Set{Symbol}()
    slotnames  = Symbol[]
    for name in framecode.code.slotnames
        if name ∉ uslotnames
            push!(slotnames, name)
            push!(uslotnames, name)
        end
    end
    assignments = Expr[]
    framename = gensym("frame")
    default = Unassigned()
    for i = 1:length(slotnames)
        slotname = framecode.code.slotnames[i]
        qslotname = QuoteNode(slotname)
        getexpr = :(something($framename.locals[$framename.last_reference[$qslotname]]))
        push!(assignments, Expr(:(=), slotname, :(haskey($framename.last_reference, $qslotname) ? $getexpr : $default)))
    end
    if ismeth
        syms = sparam_syms(framecode.scope)
        for i = 1:length(syms)
            push!(assignments, Expr(:(=), syms[i], :($framename.sparams[$i])))
        end
    end
    funcname = ismeth ? gensym("slotfunction") : gensym(Symbol(framecode.scope.name, "_slotfunction"))
    return Expr(:function, Expr(:call, funcname, framename), Expr(:block, assignments..., body))
end

const Condition = Union{Nothing,Expr,Tuple{Module,Expr}}
_unpack(condition) = isa(condition, Expr) ? (Main, condition) : condition

## The fundamental implementations of breakpoint-setting
function breakpoint!(framecode::JuliaFrameCode, pc, condition::Condition=nothing)
    stmtidx = convert(Int, pc)
    if condition === nothing
        framecode.breakpoints[stmtidx] = BreakpointState()
    else
        mod, cond = _unpack(condition)
        fex = prepare_slotfunction(framecode, cond)
        framecode.breakpoints[stmtidx] = BreakpointState(true, Core.eval(mod, fex))
    end
    return add_breakpoint(framecode, stmtidx)
end
breakpoint!(frame::JuliaStackFrame, pc=frame.pc[], condition::Condition=nothing) =
    breakpoint!(frame.code, pc, condition)

"""
    enable(bp::BreakpointRef)

Enable breakpoint `bp`.
"""
enable(bp::BreakpointRef)  = bp[] = true

"""
    disable(bp::BreakpointRef)

Disable breakpoint `bp`. Disabled breakpoints can be re-enabled with [`enable`](@ref).
"""
disable(bp::BreakpointRef) = bp[] = false

"""
    remove(bp::BreakpointRef)

Remove (delete) breakpoint `bp`. Removed breakpoints cannot be re-enabled.
"""
function remove(bp::BreakpointRef)
    idx = findfirst(isequal(bp), _breakpoints)
    deleteat!(_breakpoints, idx)
    bp.framecode.breakpoints[bp.stmtidx] = BreakpointState(false, falsecondition)
    return nothing
end

"""
    enable()

Enable all breakpoints.
"""
enable() = for bp in _breakpoints enable(bp) end

"""
    disable()

Disable all breakpoints.
"""
disable() = for bp in _breakpoints disable(bp) end

"""
    remove()

Remove all breakpoints.
"""
function remove()
    for bp in _breakpoints
        bp.framecode.breakpoints[bp.stmtidx] = BreakpointState(false, falsecondition)
    end
    empty!(_breakpoints)
    return nothing
end

"""
    breakpoint(f, sig)
    breakpoint(f, sig, line)
    breakpoint(f, sig, condition)
    breakpoint(f, sig, line, condition)
    breakpoint(...; enter_generated=false)

Add a breakpoint to `f` with the specified argument types `sig`.
Optionally specify an absolute line number `line` in the source file; the default
is to break upon entry at the first line of the body.
Without `condition`, the breakpoint will be triggered every time it is encountered;
the second only if `condition` evaluates to `true`.
`condition` should be written in terms of the arguments and local variables of `f`.

# Example
```julia
function radius2(x, y)
    return x^2 + y^2
end

breakpoint(radius2, Tuple{Int,Int}, :(y > x))
```
"""
function breakpoint(f, sig::Type, line::Integer, condition::Condition=nothing; enter_generated=false)
    method = which(f, sig)
    framecode, _ = prepare_framecode(method, sig; enter_generated=enter_generated)
    # Don't use statementnumber(method, line) in case it's enter_generated
    linec = line - whereis(method)[2] + method.line
    stmtidx = statementnumber(framecode, linec)
    breakpoint!(framecode, stmtidx, condition)
end
function breakpoint(f, sig::Type, condition::Condition=nothing; enter_generated=false)
    method = which(f, sig)
    framecode, _ = prepare_framecode(method, sig; enter_generated=enter_generated)
    breakpoint!(framecode, 1, condition)
end

"""
    breakpoint(method::Method)
    breakpoint(method::Method, line)
    breakpoint(method::Method, condition::Expr)
    breakpoint(method::Method, line, condition::Expr)

Add a breakpoint to `method`.
"""
function breakpoint(method::Method, line::Integer, condition::Condition=nothing)
    framecode, stmtidx = statementnumber(method, line)
    breakpoint!(framecode, stmtidx, condition)
end
function breakpoint(method::Method, condition::Condition=nothing)
    framecode = get_framecode(method)
    breakpoint!(framecode, 1, condition)
end

"""
    breakpoint(f)
    breakpoint(f, condition)

Break-on-entry to all methods of `f`.
"""
function breakpoint(f, condition::Condition=nothing)
    bps = BreakpointRef[]
    for method in methods(f)
        push!(bps, breakpoint(method, condition))
    end
    return bps
end

"""
    @breakpoint f(args...) condition=nothing
    @breakpoint f(args...) line condition=nothing

Break upon entry, or at the specified line number, in the method called by `f(args...)`.
Optionally supply a condition expressed in terms of the arguments and internal variables
of the method.
If `line` is supplied, it must be a literal integer.

# Example

Suppose a method `mysum` is defined as follows, where the numbers to the left are the line
number in the file:

```
12 function mysum(A)
13     s = zero(eltype(A))
14     for a in A
15         s += a
16     end
17     return s
18 end
```

Then

```
@breakpoint mysum(A) 15 s>10
```

would cause execution of the loop to break whenever `s>10`.
"""
macro breakpoint(call_expr, args...)
    whichexpr = InteractiveUtils.gen_call_with_extracted_types(__module__, :which, call_expr)
    haveline, line, condition = false, 0, nothing
    while !isempty(args)
        arg = first(args)
        if isa(arg, Integer)
            haveline, line = true, arg
        else
            condition = arg
        end
        args = Base.tail(args)
    end
    condexpr = condition === nothing ? nothing : Expr(:quote, condition)
    if haveline
        return quote
            local method = $whichexpr
            $breakpoint(method, $line, $condexpr)
        end
    else
        return quote
            local method = $whichexpr
            $breakpoint(method, $condexpr)
        end
    end
end

end


text:
module PolynomialBasis

using LinearAlgebra
import DynamicPolynomials
import StaticPolynomials

DP = DynamicPolynomials
SP = StaticPolynomials

include("lagrange_polynomials.jl")
include("hermite_polynomials.jl")
include("common_definitions.jl")
include("lagrange_basis.jl")
include("hermite_basis.jl")
include("lagrange_tensor_product_basis.jl")
include("hermite_tensor_product_basis.jl")
include("interpolation.jl")

export LagrangeBasis, CubicHermiteBasis, LagrangeTensorProductBasis,
        HermiteTensorProductBasis, gradient, hessian,
        InterpolatingPolynomial, update!, interpolation_points, basis,
        number_of_basis_functions

end # module


text:
types = [Bool, Char, Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64]

for t in types
    println("For type ", lpad(t,6), " size is $(sizeof(t)) 8-bit bytes, or ",
        lpad(string(8*sizeof(t)), 2), " bits.")
end

primitive type MyInt24 24 end

println("\nFor the 24-bit user defined type MyInt24, size is ", sizeof(MyInt24), " bytes.")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Functions for initializing HMM parameters from the observations

function kmeans_init!(hmm::AbstractHMM, observations; kwargs...)
    K = size(hmm, 1)

    res = kmeans(permutedims(observations), size(hmm, 1); kwargs...)
    seq = res.assignments

    # Initialize A
    copyto!(hmm.A, gettransmat(seq)[2])
    @check istransmat(hmm.A)

    # Initialize B
    for i in OneTo(K)
        observations_ = view(observations, seq .== i, :)
        if length(observations_) > 0
            hmm.B[i] = fit_mle(typeof(hmm.B[i]), permutedims(observations_))
        end
    end
end

max_stars_repo_name: Shu-Lea-Lai/module-MultiHMMBase


text:
"""
BigFloat operations modifying a buffer. The first argument is overwritten by the result.
All operations of MPRF returning a BigFloat are simulated.
According to the MPFR C-functions, the result variable may be identical to one of the arguments.
"""

module MPFR_MODIFY

export
    dot!, rotate!, neg!,
    mul!, add!, sub!, dv!, pow!, ceil!, copysign!, div!,
    exp!, exp2!, exponent!, factorial!, floor!, fma!, fms!, hypot!,
    ldexp!, log!, log2!, log10!, max!, min!, mod!, modf!,
    nextfloat!, prevfloat!, rem!, rem2pi!, round!,
    sqrt!, trunc!, exp10!, expm1!,
    gamma!, lgamma!, log1p!,
    sin!, cos!, tan!, sec!, csc!, cot!, acos!, asin!, atan!,
    cosh!, sinh!, tanh!, sech!, csch!, coth!, acosh!, asinh!, atanh!, atan2!,
    cbrt!, significand, frexp!

import Base: broadcast!, copy!, sum!, A_mul_B!
import Base.GMP: ClongMax, CulongMax, CdoubleMax, Limb
import Base.MPFR.ROUNDING_MODE

# Basic arithmetic without promotion
for (fJ, fC) in ((:add!,:add), (:mul!,:mul))
    @eval begin
        # BigFloat
        function ($fJ)(z::BigFloat, x::BigFloat, y::BigFloat)
            ccall(($(string(:mpfr_,fC)),:libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, &y, ROUNDING_MODE[])
            return z
        end

        # Unsigned Integer
        function ($fJ)(z::BigFloat, x::BigFloat, c::CulongMax)
            ccall(($(string(:mpfr_,fC,:_ui)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Culong, Int32), &z, &x, c, ROUNDING_MODE[])
            return z
        end
        ($fJ)(z::BigFloat, c::CulongMax, x::BigFloat) = ($fJ)(z,x,c)

        # Signed Integer
        function ($fJ)(z::BigFloat, x::BigFloat, c::ClongMax)
            ccall(($(string(:mpfr_,fC,:_si)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Clong, Int32), &z, &x, c, ROUNDING_MODE[])
            return z
        end
        ($fJ)(z::BigFloat, c::ClongMax, x::BigFloat) = ($fJ)(z,x,c)

        # Float32/Float64
        function ($fJ)(z::BigFloat, x::BigFloat, c::CdoubleMax)
            ccall(($(string(:mpfr_,fC,:_d)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Cdouble, Int32), &z, &x, c, ROUNDING_MODE[])
            return z
        end
        ($fJ)(z::BigFloat, c::CdoubleMax, x::BigFloat) = ($fJ)(z,x,c)

        # BigInt
        function ($fJ)(z::BigFloat, x::BigFloat, c::BigInt)
            ccall(($(string(:mpfr_,fC,:_z)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigInt}, Int32), &z, &x, &c, ROUNDING_MODE[])
            return z
        end
        ($fJ)(z::BigFloat, c::BigInt, x::BigFloat) = ($fJ)(z,x,c)
    end
end

for (fJ, fC) in ((:sub!,:sub), (:dv!,:div))
    @eval begin
        # BigFloat
        function ($fJ)(z::BigFloat, x::BigFloat, y::BigFloat)
            ccall(($(string(:mpfr_,fC)),:libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, &y, ROUNDING_MODE[])
            return z
        end

        # Unsigned Int
        function ($fJ)(z::BigFloat, x::BigFloat, c::CulongMax)
            ccall(($(string(:mpfr_,fC,:_ui)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Culong, Int32), &z, &x, c, ROUNDING_MODE[])
            return z
        end
        function ($fJ)(z::BigFloat, c::CulongMax, x::BigFloat)
            ccall(($(string(:mpfr_,:ui_,fC)), :libmpfr), Int32, (Ptr{BigFloat}, Culong, Ptr{BigFloat}, Int32), &z, c, &x, ROUNDING_MODE[])
            return z
        end

        # Signed Integer
        function ($fJ)(z::BigFloat, x::BigFloat, c::ClongMax)
            ccall(($(string(:mpfr_,fC,:_si)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Clong, Int32), &z, &x, c, ROUNDING_MODE[])
            return z
        end
        function ($fJ)(z::BigFloat, c::ClongMax, x::BigFloat)
            ccall(($(string(:mpfr_,:si_,fC)), :libmpfr), Int32, (Ptr{BigFloat}, Clong, Ptr{BigFloat}, Int32), &z, c, &x, ROUNDING_MODE[])
            return z
        end

        # Float32/Float64
        function ($fJ)(z::BigFloat, x::BigFloat, c::CdoubleMax)
            ccall(($(string(:mpfr_,fC,:_d)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Cdouble, Int32), &z, &x, c, ROUNDING_MODE[])
            return z
        end
        function ($fJ)(z::BigFloat, c::CdoubleMax, x::BigFloat)
            ccall(($(string(:mpfr_,:d_,fC)), :libmpfr), Int32, (Ptr{BigFloat}, Cdouble, Ptr{BigFloat}, Int32), &z, c, &x, ROUNDING_MODE[])
            return z
        end

        # BigInt
        function ($fJ)(z::BigFloat, x::BigFloat, c::BigInt)
            ccall(($(string(:mpfr_,fC,:_z)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigInt}, Int32), &z, &x, &c, ROUNDING_MODE[])
            return z
        end
        # no :mpfr_z_div function
    end
end

function sub!(z::BigFloat, c::BigInt, x::BigFloat)
    ccall((:mpfr_z_sub, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigInt}, Ptr{BigFloat}, Int32), &z, &c, &x, ROUNDING_MODE[])
    return z
end

function fma!(r::BigFloat, x::BigFloat, y::BigFloat, z::BigFloat)
    ccall(("mpfr_fma",:libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &r, &x, &y, &z, ROUNDING_MODE[])
    return r
end

function fms!(r::BigFloat, x::BigFloat, y::BigFloat, z::BigFloat)
    ccall(("mpfr_fms",:libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &r, &x, &y, &z, ROUNDING_MODE[])
    return r
end

# div
# BigFloat
function div!(z::BigFloat, x::BigFloat, y::BigFloat)
    ccall((:mpfr_div,:libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, &y, to_mpfr(RoundToZero))
    ccall((:mpfr_trunc, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}), &z, &z)
    return z
end

# Unsigned Int
function div!(z::BigFloat, x::BigFloat, c::CulongMax)
    ccall((:mpfr_div_ui, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Culong, Int32), &z, &x, c, to_mpfr(RoundToZero))
    ccall((:mpfr_trunc, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}), &z, &z)
    return z
end
function div!(z::BigFloat, c::CulongMax, x::BigFloat)
    ccall((:mpfr_ui_div, :libmpfr), Int32, (Ptr{BigFloat}, Culong, Ptr{BigFloat}, Int32), &z, c, &x, to_mpfr(RoundToZero))
    ccall((:mpfr_trunc, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}), &z, &z)
    return z
end

# Signed Integer
function div!(z::BigFloat, x::BigFloat, c::ClongMax)
    ccall((:mpfr_div_si, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Clong, Int32), &z, &x, c, to_mpfr(RoundToZero))
    ccall((:mpfr_trunc, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}), &z, &z)
    return z
end
function div!(z::BigFloat, c::ClongMax, x::BigFloat)
    ccall((:mpfr_si_div, :libmpfr), Int32, (Ptr{BigFloat}, Clong, Ptr{BigFloat}, Int32), &z, c, &x, to_mpfr(RoundToZero))
    ccall((:mpfr_trunc, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}), &z, &z)
    return z
end

# Float32/Float64
function div!(z::BigFloat, x::BigFloat, c::CdoubleMax)
    ccall((:mpfr_div_d, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Cdouble, Int32), &z, &x, c, to_mpfr(RoundToZero))
    ccall((:mpfr_trunc, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}), &z, &z)
    return z
end
function div!(z::BigFloat, c::CdoubleMax, x::BigFloat)
    ccall((:mpfr_d_div, :libmpfr), Int32, (Ptr{BigFloat}, Cdouble, Ptr{BigFloat}, Int32), &z, c, &x, to_mpfr(RoundToZero))
    ccall((:mpfr_trunc, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}), &z, &z)
    return z
end

# BigInt
function div!(z::BigFloat, x::BigFloat, c::BigInt)
    ccall((:mpfr_div_z, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigInt}, Int32), &z, &x, &c, to_mpfr(RoundToZero))
    ccall((:mpfr_trunc, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}), &z, &z)
    return z
end

# More efficient commutative operations
for (fJ, fC, fI) in ((:add!, :add, 0), (:mul!, :mul, 1))
    @eval begin
        function ($fJ)(z::BigFloat, a::BigFloat, b::BigFloat, c::BigFloat)
            ccall(($(string(:mpfr_,fC)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &a, &b, ROUNDING_MODE[])
            ccall(($(string(:mpfr_,fC)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &z, &c, ROUNDING_MODE[])
            return z
        end
        function ($fJ)(z::BigFloat, a::BigFloat, b::BigFloat, c::BigFloat, d::BigFloat)
            ccall(($(string(:mpfr_,fC)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &a, &b, ROUNDING_MODE[])
            ccall(($(string(:mpfr_,fC)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &z, &c, ROUNDING_MODE[])
            ccall(($(string(:mpfr_,fC)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &z, &d, ROUNDING_MODE[])
            return z
        end
        function ($fJ)(z::BigFloat, a::BigFloat, b::BigFloat, c::BigFloat, d::BigFloat, e::BigFloat)
            ccall(($(string(:mpfr_,fC)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &a, &b, ROUNDING_MODE[])
            ccall(($(string(:mpfr_,fC)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &z, &c, ROUNDING_MODE[])
            ccall(($(string(:mpfr_,fC)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &z, &d, ROUNDING_MODE[])
            ccall(($(string(:mpfr_,fC)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &z, &e, ROUNDING_MODE[])
            return z
        end
    end
end

function neg!(z::BigFloat)
    ccall((:mpfr_neg, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &z, ROUNDING_MODE[])
    return z
end

function neg!(z::BigFloat, x::BigFloat)
    ccall((:mpfr_neg, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[])
    return z
end

function sqrt!(z::BigFloat, x::BigFloat)
    isnan(x) && return x
    ccall((:mpfr_sqrt, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[])
    if isnan(z)
        throw(DomainError())
    end
    return z
end

sqrt!(z::BigFloat, x::BigInt) = sqrt!(z, BigFloat(x))

function pow!(z::BigFloat, x::BigFloat, y::BigFloat)
    ccall((:mpfr_pow, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, &y, ROUNDING_MODE[])
    return z
end

function pow!(z::BigFloat, x::BigFloat, y::CulongMax)
    ccall((:mpfr_pow_ui, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Culong, Int32), &z, &x, y, ROUNDING_MODE[])
    return z
end

function pow!(z::BigFloat, x::BigFloat, y::ClongMax)
    ccall((:mpfr_pow_si, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Clong, Int32), &z, &x, y, ROUNDING_MODE[])
    return z
end

function pow!(z::BigFloat, x::BigFloat, y::BigInt)
    ccall((:mpfr_pow_z, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigInt}, Int32), &z, &x, &y, ROUNDING_MODE[])
    return z
end

pow!(z::BigFloat, x::BigFloat, y::Integer)  = typemin(Clong)  <= y <= typemax(Clong)  ? pow!(z, x, Clong(y))  : pow!(z, x, BigInt(y))

pow!(z::BigFloat, x::BigFloat, y::Unsigned) = typemin(Culong) <= y <= typemax(Culong) ? pow!(z, x, Culong(y)) : pow!(z, x, BigInt(y))

for f in (:exp, :exp2, :exp10, :expm1, :cosh, :sinh, :tanh, :sech, :csch, :coth, :cbrt)
@eval function $(Symbol(f, "!"))(z::BigFloat, x::BigFloat)
        ccall(($(string(:mpfr_,f)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[])
        return z
    end
end

# return log(2)
function big_ln2!(c::BigFloat)
    ccall((:mpfr_const_log2, :libmpfr), Cint, (Ptr{BigFloat}, Int32),
          &c, MPFR.ROUNDING_MODE[])
    return c
end

function ldexp!(z::BigFloat, x::BigFloat, n::Clong)
    ccall((:mpfr_mul_2si, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Clong, Int32), &z, &x, n, ROUNDING_MODE[])
    return z
end
function ldexp!(z::BigFloat, x::BigFloat, n::Culong)
    ccall((:mpfr_mul_2ui, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Culong, Int32), &z, &x, n, ROUNDING_MODE[])
    return z
end
ldexp!(z::BigFloat, x::BigFloat, n::ClongMax) = ldexp!(z, x, convert(Clong, n))
ldexp!(z::BigFloat, x::BigFloat, n::CulongMax) = ldexp!(z, x, convert(Culong, n))
ldexp!(z::BigFloat, x::BigFloat, n::Integer) = mul!(z, x, exp2!(BigFloat(n)))

function factorial!(z::BigFloat, x::BigFloat)
    if x < 0 || !isinteger(x)
        throw(DomainError())
    end
    ui = convert(Culong, x)
    ccall((:mpfr_fac_ui, :libmpfr), Int32, (Ptr{BigFloat}, Culong, Int32), &z, ui, ROUNDING_MODE[])
    return z
end

function hypot!(z::BigFloat, x::BigFloat, y::BigFloat)
    ccall((:mpfr_hypot, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, &y, ROUNDING_MODE[])
    return z
end

for f in (:log, :log2, :log10)
  @eval function $(Symbol(f, "!"))(z::BigFloat, x::BigFloat)
        if x < 0
            throw(DomainError())
        end
        ccall(($(string(:mpfr_,f)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[])
        return z
    end
end

function log1p!(z::BigFloat, x::BigFloat)
    if x < -1
        throw(DomainError())
    end
    ccall((:mpfr_log1p, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[])
    return z
end

function max!(z::BigFloat, x::BigFloat, y::BigFloat)
    isnan(x) && return x
    isnan(y) && return y
    ccall((:mpfr_max, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, &y, ROUNDING_MODE[])
    return z
end

function min!(z::BigFloat, x::BigFloat, y::BigFloat)
    isnan(x) && return x
    isnan(y) && return y
    ccall((:mpfr_min, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, &y, ROUNDING_MODE[])
    return z
end

function modf!(zint::BigFloat, zfloat::BigFloat, x::BigFloat)
  if isinf(x)
    copy!(zint, x)
    copy!(zfloat, BigFloat(NaN))
  end
  ccall((:mpfr_modf, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &zint, &zfloat, &x, ROUNDING_MODE[])
end

function rem!(z::BigFloat, x::BigFloat, y::BigFloat)
    ccall((:mpfr_fmod, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, &y, ROUNDING_MODE[])
    return z
end

function rem!(z::BigFloat, x::BigFloat, y::BigFloat, ::RoundingMode{:Nearest})
    ccall((:mpfr_remainder, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, &y, ROUNDING_MODE[])
    return z
end

# TODO: use a higher-precision BigFloat(pi) here?
function rem2pi!(z::BigFloat, x::BigFloat, r::RoundingMode)
  mul!(z, BigFloat(pi), 2)
  rem!(z, x, z, r)
end

function sum!(z::BigFloat, arr::AbstractArray{BigFloat})
copy!(z, BigFloat(0))
    for i in arr
        ccall((:mpfr_add, :libmpfr), Int32,
            (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Cint),
            &z, &z, &i, 0)
    end
    return z
end

# Functions for which NaN results are converted to DomainError, following Base
for f in (:sin,:cos,:tan,:sec,:csc,
          :acos,:asin,:atan,:acosh,:asinh,:atanh, :gamma)
  @eval begin
    function $(Symbol(f, "!"))(z::BigFloat, x::BigFloat)
            if isnan(x)
                return x
            end
            ccall(($(string(:mpfr_,f)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[])
            if isnan(z)
                throw(DomainError())
            end
            return z
        end
    end
end

# log of absolute value of gamma function
const lgamma_signp2 = Ref{Cint}()
function lgamma!(z::BigFloat, x::BigFloat)
    ccall((:mpfr_lgamma,:libmpfr), Cint, (Ptr{BigFloat}, Ptr{Cint}, Ptr{BigFloat}, Int32), &z, lgamma_signp2, &x, ROUNDING_MODE[])
    return z
end

lgamma_r!(z::BigFloat, x::BigFloat) = (lgamma!(z, x), lgamma_signp2[])

function atan2!(z::BigFloat, y::BigFloat, x::BigFloat)
    ccall((:mpfr_atan2, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &y, &x, ROUNDING_MODE[])
    return z
end

function copysign!(z::BigFloat, x::BigFloat, y::BigFloat)
    ccall((:mpfr_copysign, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, &y, ROUNDING_MODE[])
    return z
end

function frexp!(z::BigFloat, x::BigFloat)
    c = Ref{Clong}()
    ccall((:mpfr_frexp, :libmpfr), Int32, (Ptr{Clong}, Ptr{BigFloat}, Ptr{BigFloat}, Cint), c, &z, &x, ROUNDING_MODE[])
    return (z, c[])
end

function significand!(z::BigFloat, x::BigFloat)
    c = Ref{Clong}()
    ccall((:mpfr_frexp, :libmpfr), Int32, (Ptr{Clong}, Ptr{BigFloat}, Ptr{BigFloat}, Cint), c, &z, &x, ROUNDING_MODE[])
    # Double the significand to make it work as Base.significand
    ccall((:mpfr_mul_si, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Clong, Int32), &z, &z, 2, ROUNDING_MODE[])
    return z
end

function round!(z::BigFloat, x::BigFloat)
    ccall((:mpfr_rint, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Cint), &z, &x, ROUNDING_MODE[])
    return z
end
function round!(z::BigFloat, x::BigFloat,::RoundingMode{:NearestTiesAway})
    ccall((:mpfr_round, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}), &z, &x)
    return z
end

function copy!(z::BigFloat, x::BigFloat)
    ccall((:mpfr_set, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32),
          &z, &x, ROUNDING_MODE[])
    return z
end

function copy!(z::BigFloat, x::ClongMax)
    ccall((:mpfr_set_si, :libmpfr), Int32, (Ptr{BigFloat}, Clong, Int32),
          &z, x, ROUNDING_MODE[])
    return z
end

function copy!(z::BigFloat, x::CulongMax)
    ccall((:mpfr_set_ui, :libmpfr), Int32, (Ptr{BigFloat}, Culong, Int32),
          &z, x, ROUNDING_MODE[])
    return z
end

function copy!(z::BigFloat, x::CdoubleMax)
    ccall((:mpfr_set_d, :libmpfr), Int32, (Ptr{BigFloat}, Cdouble, Int32),
          &z, x, ROUNDING_MODE[])
    return z
end

function copy!(z::BigFloat, x::BigInt)
    ccall((:mpfr_set_z, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigInt}, Int32),
          &z, &x, ROUNDING_MODE[])
    return z
end

copy!(z::BigFloat, x::Real) = copy!(z, BigFloat(x))

function nextfloat!(z::BigFloat, x::BigFloat)
  copy!(z, x)
  ccall((:mpfr_nextabove, :libmpfr), Int32, (Ptr{BigFloat},), &z) != 0
  return z
end

function prevfloat!(z::BigFloat, x::BigFloat)
  copy!(z, x)
  ccall((:mpfr_nextbelow, :libmpfr), Int32, (Ptr{BigFloat},), &z) != 0
  return z
end

# broadcast!(::typeof(identity), z::BigFloat, x::Real) = copy!(z, x)

function checklength(n::Int, m::Int)
  n != m &&
    throw(DimensionMismatch("first array has length $(n) which does not match the length of the second, $(m)."))
  n
end

checklength(a, b) = checklength(length(a), length(b))
checkmult(a::AbstractArray{T,2}, b) where T = checklength(size(a, 2), length(b))
checkmult(a, b::AbstractArray{T,2}) where T = checklenght(length(a), size(b, 2))

function dot!(z::BigFloat, a::AbstractArray{T,1}, b::AbstractArray{U,1}) where T <: Real where U <: Real

  n = checklength(a, b)
  _dot!(z, n, a, b)
end

function _dot!(z::BigFloat, n::Int, a::AbstractArray{T,1}, b::AbstractArray{T,1}) where T <: BigFloat
  copy!(z, 0)
  for i in 1:n
    fma!(z, a[i], b[i], z)
  end
  z
end

function _dot!(z::BigFloat, n::Int, a::AbstractArray{T,1}, b::AbstractArray{U,1}) where T <: BigFloat where U <: Real

  copy!(z, 0)
  c = BigFloat()
  for i in 1:n
    copy!(c, b[i])
    fma!(z, a[i], c, z)
  end
  z
end

function _dot!(z::BigFloat, n::Int, a::AbstractArray{T,1}, b::AbstractArray{U,1}) where T <: Real where U <: BigFloat
  _dot!(z, n, b, a)
end

function _dot!(z::BigFloat, n::Int, a::AbstractArray{T,1}, b::AbstractArray{U,1}) where T <: Real where U <: Real

  copy!(z, 0)
  c = BigFloat()
  for i in 1:n
    copy!(c, a[i])
    mul!(c, c, b[i])
    add!(z, z, c)
  end
  z
end

function mul!(z::AbstractVector{BigFloat}, a::AbstractArray{BigFloat,2}, b::AbstractVector{BigFloat})

  checkmult(a, b)
  n = checkmult(z, a)
  for i = 1:n
    isassigned(z, i) || ( z[i] = BigFloat() )
    dot!(z[i], view(a, i,:), b)
  end
  z
end

function mul!(z::AbstractVector{BigFloat}, a::AbstractVector{BigFloat}, b::Real)
  n = checklength(z, a)
  for i = 1:n
    isassigned(z, i) || ( z[i] = BigFloat() )
    mul!(z[i], a[i], b)
  end
  z
end

function add!(z::AbstractVector{BigFloat}, a::AbstractVector{BigFloat}, b::AbstractVector{<:Real})
  n = checklength(z, a)
  for i = 1:n
    isassigned(z, i) || ( z[i] = BigFloat() )
    add!(z[i], a[i], b[i])
  end
  z
end

function sub!(z::AbstractVector{BigFloat}, a::AbstractVector{BigFloat}, b::AbstractVector{<:Real})
  n = checklength(z, a)
  for i = 1:n
    isassigned(z, i) || ( z[i] = BigFloat() )
    sub!(z[i], a[i], b[i])
  end
  z
end

"""
Apply givens rotation, defined by `(c, s)` to `[a1, a2]`.
Use temporary work space `t[1:2]`. No aliasing of any used varables allowed.
"""
function rotate!(a1::BigFloat, a2::BigFloat, c::BigFloat, s::BigFloat, t::AbstractVector{BigFloat})

  mul!(t[1], a1, c)
  mul!(t[2], a1, s)
  fma!(a1, a2, s, t[1])
  fms!(a2, a2, c, t[2])
  nothing
end

end # module

max_stars_repo_name: KlausC/LinAlgBigFloat

max_issues_repo_name: KlausC/LinAlgBigFloat

max_forks_repo_name: KlausC/LinAlgBigFloat


text:
module CoreEcoliModel

# internal includes -
include("Include.jl")

#
# maximize_specific_growth_rate(path_to_measurements_file::String; number_of_samples::Int64 = 100)
export maximize_specific_growth_rate
export maximize_flux_at_index
export moma_calculation
export flux_variabilty_analysis

# Data dictionary methods -
export generate_default_data_dictionary
export constrain_measured_fluxes
export constrain_measured_metabolites
export optimize_specific_growth_rate

# Export utility methods -
export write_model_reactions_to_file

# export types -
export VLOptimalFluxResult

end # module


text:
using Ipopt
using MathOptInterface
const MOI = MathOptInterface

#Boilerplate setup code to interface with IPOPT.

struct ProblemMOI <: MOI.AbstractNLPEvaluator
    n_nlp::Int
    m_nlp::Int
    idx_ineq
    obj_grad::Bool
    con_jac::Bool
    sparsity_jac
    sparsity_hess
    primal_bounds
    constraint_bounds
    hessian_lagrangian::Bool
end

function ProblemMOI(n_nlp,m_nlp;
        idx_ineq=(1:0),
        obj_grad=true,
        con_jac=true,
        sparsity_jac=sparsity_jacobian(n_nlp,m_nlp),
        sparsity_hess=sparsity_hessian(n_nlp,m_nlp),
        primal_bounds=primal_bounds(n_nlp),
        constraint_bounds=constraint_bounds(m_nlp,idx_ineq=idx_ineq),
        hessian_lagrangian=false)

    ProblemMOI(n_nlp,m_nlp,
        idx_ineq,
        obj_grad,
        con_jac,
        sparsity_jac,
        sparsity_hess,
        primal_bounds,
        constraint_bounds,
        hessian_lagrangian)
end

function constraint_bounds(m; idx_ineq=(1:0))
    c_l = zeros(m)

    c_u = zeros(m)
    c_u[idx_ineq] .= Inf

    return c_l, c_u
end

function row_col!(row,col,r,c)
    for cc in c
        for rr in r
            push!(row,convert(Int,rr))
            push!(col,convert(Int,cc))
        end
    end
    return row, col
end

function sparsity_jacobian(n,m)

    row = []
    col = []

    r = 1:m
    c = 1:n

    row_col!(row,col,r,c)

    return collect(zip(row,col))
end

function sparsity_hessian(n,m)

    row = []
    col = []

    r = 1:m
    c = 1:n

    row_col!(row,col,r,c)

    return collect(zip(row,col))
end

function MOI.eval_objective(prob::MOI.AbstractNLPEvaluator, x)
    objective(x)
end

function MOI.eval_objective_gradient(prob::MOI.AbstractNLPEvaluator, grad_f, x)
    ForwardDiff.gradient!(grad_f,objective,x)
    return nothing
end

function MOI.eval_constraint(prob::MOI.AbstractNLPEvaluator,g,x)
    constraint!(g,x)
    return nothing
end

function MOI.eval_constraint_jacobian(prob::MOI.AbstractNLPEvaluator, jac, x)
    ForwardDiff.jacobian!(reshape(jac,prob.m_nlp,prob.n_nlp), constraint!, zeros(prob.m_nlp), x)
    return nothing
end

function MOI.features_available(prob::MOI.AbstractNLPEvaluator)
    return [:Grad, :Jac]
end

MOI.initialize(prob::MOI.AbstractNLPEvaluator, features) = nothing
MOI.jacobian_structure(prob::MOI.AbstractNLPEvaluator) = prob.sparsity_jac

function ipopt_solve(x0,prob::MOI.AbstractNLPEvaluator;
    tol=1.0e-3,c_tol=1.0e-3,max_iter=10000, print=0)
    x_l, x_u = prob.primal_bounds
    c_l, c_u = prob.constraint_bounds

    nlp_bounds = MOI.NLPBoundsPair.(c_l,c_u)
    block_data = MOI.NLPBlockData(nlp_bounds,prob,true)

    solver = Ipopt.Optimizer()
    solver.options["max_iter"] = max_iter
    solver.options["tol"] = tol
    solver.options["constr_viol_tol"] = c_tol

    #Uncomment the following line to turn off verbose IPOPT output
    solver.options["print_level"] = print

    x = MOI.add_variables(solver,prob.n_nlp)

    for i = 1:prob.n_nlp
        xi = MOI.SingleVariable(x[i])
        # xi = MOI.x[i]
        MOI.add_constraint(solver, xi, MOI.LessThan(x_u[i]))
        MOI.add_constraint(solver, xi, MOI.GreaterThan(x_l[i]))
        MOI.set(solver, MOI.VariablePrimalStart(), x[i], x0[i])
    end

    # Solve the problem
    MOI.set(solver, MOI.NLPBlock(), block_data)
    MOI.set(solver, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.optimize!(solver)

    # Get the solution
    res = MOI.get(solver, MOI.VariablePrimal(), x)

    return res
end

max_stars_repo_name: bbokser/hopper-sim

max_issues_repo_name: bbokser/hopper-sim

max_forks_repo_name: bbokser/hopper-sim


text:
using CompScienceMeshes, BEAST

o, x, y, z = euclidianbasis(3)

sol = 5.0;
Δt, Nt = 100.0/sol,200

D, Δx = 1.0, 0.45
Γ = meshsphere(D, Δx)
X = raviartthomas(Γ)

(A, b, c) = butcher_tableau_radau_2stages();
V = StagedTimeStep(X, c, Δt, Nt);

duration, delay, amplitude = 2000.0/sol, 2000.0/sol, 1.0
gaussian = creategaussian(duration, delay, duration)

direction, polarisation = z, x
E = planewave(polarisation, direction, derive(gaussian), sol)

LaplaceEFIO(s::T) where {T} = MWSingleLayer3D(-s/sol, s*s/sol, T(sol));
kmax = 15;
rho = 1.0001;
T = RungeKuttaConvolutionQuadrature(LaplaceEFIO, A, b, Δt, kmax, rho);

@hilbertspace j
@hilbertspace j′
tdefie = @discretise T[j′,j] == -1E[j′]   j∈V  j′∈V
xefie_irk = solve(tdefie)


text:
function GRASP(cost, matrix, n, m, alpha)
#    println("%%%%%%% Debut GRASP %%%%%%%")
#    println(alpha)

    #Creating a set of lines that will be avaluated
    desactive_condition= zeros(Bool, m) #if desactive_condition[j]=0 the line will be avaluated
    stop1= ones(Bool,m)

    #Creation the an array of activated variables
    variables_actives = ones(Bool, n)
    stop2= zeros(Bool , n)
#    println("\n variables_actives: ", variables_actives, "\n stop2: ", stop2)

    #creating my utility array for my m variables
    util = zeros(Float64, n)

    #initialisation de la Solution
    SOL = zeros(Float64, n)

#    println( "\n Matrix: ", matrix, "\n Cost: ", cost,"\n n= ", n, "   m= ", m,"\n desactive_condition: ", desactive_condition, "\n stop: ", stop1, "\n util: ", util, "\n SOL: ", SOL, "\n \n")

#    ite = 0
    #Création de la solution
    while desactive_condition!=stop1 && variables_actives!=stop2

#      println("\n ++++++++++++++++++++++++++++++++++++++++")
#      println("\n itération: ", ite)
      #println("\n desactive_condition : ", desactive_condition, "\n stop1 : " , "\n variables_actives : " , variables_actives , "\n stop2 : ")
#      println("\n ++++++++++++++++++++++++++++++++++++++++")

        #Foction d'utilité

        util = Utilite(cost, matrix, desactive_condition, n, m , variables_actives, stop1, stop2)

        #Choix parmis les candidats de util
#        PosCandidat::Int64
#        PosCandidat = PosMax(util,n)
        PosCandidat::Int64 = choixcandidat(util, alpha)

        #On ajoute le candidat à la solution
        SOL[PosCandidat] = 1

        #Desactive! le candidat selectionné
        Desactive!(PosCandidat, matrix, desactive_condition, m, variables_actives,n)
#        ite=ite+1
    end
    #println(calculz(SOL,cost,n))
    #return(SOL)
    Z = calculz(SOL,cost,n)
#    println("Solution GRASP: " , Z)
#    return(SOL, desactive_condition, Z)
    return(SOL, Z)

end

# =========================================================================== #

#Detremine an utility  based on an active (desactive_condition) set of matrix's lines
function Utilite(cost, matrix, desactive_condition, n, m , variables_actives, stop1, stop2)
#    println("==========================")
    util = zeros(Float64, n) #réinitialisation du vecteur

    #On each column of matrix
    for j=1:n
            if variables_actives[j] == 1
                    #For each line
                    for i=1:m
                        #checking if the line is active
                        if desactive_condition[i]==0 #&& i <= m
                            K=matrix[i,j]
                            util[j]=util[j]+K
#                            println("====== \n", "i: ", i ,"    j: ", j , "   k: ", K )
#                            println("\n variables_actives: ", variables_actives)
#                            println("\n desactive_condition : ", desactive_condition)
#                            println("\n util : ", util)
                        end
                    end

                # dividing the number if iterations of the variable j by its cost (not a zero)

                if util[j]!=false #On évide de divisé par 0...
#                    println("\n util avant modif: ", util)
                    util[j] =  cost[j] / util[j]
#                    println("\n util après modif: ", util)

                else
                    variables_actives[j] = 0
                end
            end
    end
#    println("\n Util: ", util)
#    println("\n \n ===================Utilite fini =============")
    return util
end

# =========================================================================== #
#=
function PosMax(util,n) #On choisie un candidat parmie les différents coûts, par ordre lexicographique
#println("%%%%%%%%%%%%%%%%%%%%")
#    Pos::Int64
    Pos=0
    ValCan=0.0

    for i=1:n

        if util[i] > ValCan

            ValCan = util[i]
            Pos = i
        end
    end
    #println("\n \n %%%%%% Pos max fini %%%%% \n pos: ", Pos)
    return Pos
end
=#

# =========================================================================== #

function Desactive!(PosCandidat, matrix, desactive_condition, m , variables_actives , n) #On désactive les lignes où est le candidat
#println("\n======================= desactive ==============================")
#println("~~~~~~~~~~~~~~~~~~~~~~~~")
    for i=1:m
#        println("\n Boucle i: ", i)
        if matrix[i,PosCandidat] == 1
#            println("\n If matrix[i,PosCandidat] == 1: " , matrix[i,PosCandidat])
            desactive_condition[i] = 1
#            println("\n desactive_condition[i] = 1: ", desactive_condition[i])
            for j=1:n
#                println("\n for j=1:n ; j: ", j)
                if matrix[i,j] ==1
#                    println("\n if matrix[i,j] ==1 :", matrix[i,j])
                    variables_actives[j] = 0
#                    println("\n variables_actives[j] = 0: ", variables_actives)
                end
            end
        end
    end
#    println("\n \n ~~~~~~~~~Desactive fini! ~~~~~~~\n ")
#    println("variables_actives: ", variables_actives)
#    println("desactive_condition: ", desactive_condition)
end

#=========================================================================#

function calculz(x,costs,m)
    z = 0
    for i in 1:m
        z+= x[i]*costs[i]
    end
    return z
end

#==========================================#
function choixcandidat(util, alpha)
    max = maximum(util)
    min = minimum(util)
#    random_utile = Array{Float64}[]
    random_utile = Float64[]
#    println("util: ", util, "\n type de util: ", typeof(util))
#    println("min: ", min, " max: ", max, " random_utile: ", random_utile)
    cpt::Int=1
    for i in 1:length(util)
#        println("debut for: ", i)
#        println("taille util: ", length(util), "   alpha: ", alpha)
        if util[i]>= min+(alpha*(max-min))
            push!(random_utile,i)
        end
    end
#    println("taille: ", length(random_utile))#, "   max: ", max, "min+(alpha*(max-min): ", min+(alpha*(max-min)))
#=    if length(random_utile)==0
        println("   max: ", max, "   min: ", min , "  min+(alpha*(max-min): ", min+(alpha*(max-min)) , "  alpha :", alpha)
        println("util: ", util)
    end
=#
    prob = rand(1:length(random_utile))
#    println("fini, random_utile: ", random_utile, " Prob: ", prob)
    return(random_utile[prob])
end

max_stars_repo_path: Archives/before the merge/DM2_1.jl

max_issues_repo_path: Archives/before the merge/DM2_1.jl

max_forks_repo_path: Archives/before the merge/DM2_1.jl


text:
##
##  Java Information Dynamics Toolkit (JIDT)
##  Copyright (C) 2012, Joseph T. Lizier
##
##  This program is free software: you can redistribute it and/or modify
##  it under the terms of the GNU General Public License as published by
##  the Free Software Foundation, either version 3 of the License, or
##  (at your option) any later version.
##
##  This program is distributed in the hope that it will be useful,
##  but WITHOUT ANY WARRANTY; without even the implied warranty of
##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##  GNU General Public License for more details.
##
##  You should have received a copy of the GNU General Public License
##  along with this program.  If not, see <http://www.gnu.org/licenses/>.
##

# = Example 1 - Transfer entropy on binary data =

# Simple transfer entropy (TE) calculation on binary data using the discrete TE calculator:

# Import the JavaCall package:
using JavaCall;

# Change location of jar to match yours:
jarLocation = "../../infodynamics.jar";
# Start the JVM supplying classpath and heap size
#  (increase memory here if you get crashes due to not enough space)
JavaCall.init(["-Djava.class.path=$(jarLocation)", "-Xmx128M"]);

# Generate some random binary data.
sourceArray=rand(0:1, 100);
destArray = [0; sourceArray[1:99]];
sourceArray2=rand(0:1, 100);

# Create a TE calculator and run it:
teClass = @jimport infodynamics.measures.discrete.TransferEntropyCalculatorDiscrete;
teCalc=teClass((jint, jint), 2, 1);
jcall(teCalc, "initialise", Void, ()); # This is how to indicate a void return type and arguments

# We can pass simple arrays of ints directly in:
jcall(teCalc, "addObservations", Void, (Array{jint,1}, Array{jint,1}),
				     sourceArray, destArray);
result = jcall(teCalc, "computeAverageLocalOfObservations", jdouble, ());
@printf("For copied source, result should be close to 1 bit : %.4f\n", result);

jcall(teCalc, "initialise", Void, ());
jcall(teCalc, "addObservations", Void, (Array{jint,1}, Array{jint,1}),
				     sourceArray2, destArray);
result2 = jcall(teCalc, "computeAverageLocalOfObservations", jdouble, ());
@printf("For random source, result should be close to 0 bits: %.4f\n", result2);

max_issues_repo_name: jmmccracken/ECA


text:
using ExcelFiles
using Mimi
using Random
using Statistics
using Test

include("../../wip/load_empirical_dist.jl")

# function load_vector(path, range, header=false)
#     tups = collect(load(path, range, header=header))
#     name = fieldnames(tups[1])[1]   # field name of first item in NamedTuple
#     map(obj -> getfield(obj, name), tups)
# end

# function load_empirical_dist(path::AbstractString,
#                              values_range::AbstractString,
#                              probs_range::AbstractString="")
#     println("Reading from '$path'")
#     values = load_vector(path, values_range)
#     probs = probs_range == "" ? nothing : load_vector(path, probs_range)
#     d = Mimi.EmpiricalDistribution(values, probs)
#     println("returning distribution $(typeof(d))")
#     return d
# end

filename = joinpath(@__DIR__, "RB-ECS-distribution.xls")
d = load_empirical_dist(filename, "Sheet1!A2:A1001", "Sheet1!B2:B1001")

# Set the seed to get repeatable results (with some caveats...)
Random.seed!(1234567)

# Get the statistical outline of the distribution
q = quantile.(Ref(d), [0.01, 0.10, 0.25, 0.50, 0.75, 0.90, 0.99])

expected = [1.47, 1.92, 2.33, 3.01, 4.16, 5.85, 9.16]
# @info "quantiles: $q\n expected: $expected"

@test isapprox(q, expected, atol=0.01)

# Test that the EmpiricalDistribution gets saved as a SampleStore and values
# get re-used across multiple scenarios

num_scenarios = 4
num_trials = 5
output_dir = "./out/"

results = zeros(num_scenarios, num_trials)

_values = collect(1:10)
_probs = 0.1 * ones(10)

mcs_test = @defmcs begin
    p = EmpiricalDistribution(_values, _probs)
end

@defcomp test begin
    p = Parameter(default = 5)
    function run_timestep(p, v, d, t) end
end

scenario_args = [
    :num => collect(1:num_scenarios)
]

function scenario_func(mcs::MonteCarloSimulation, tup::Tuple)
    nothing
end

function post_trial_func(mcs::MonteCarloSimulation, trialnum::Int, ntimesteps::Int, tup::Tuple)
    m, = mcs.models
    scenario_num, = tup
    results[scenario_num, trialnum] = m[:test, :p]
end

m = Model()
set_dimension!(m, :time, 2000:2001)
add_comp!(m, test)
set_model!(mcs_test, m)

generate_trials!(mcs_test, num_trials; sampling = RANDOM)

run_mcs(mcs_test, num_trials;
    output_dir = output_dir,
    scenario_args = scenario_args,
    scenario_func = scenario_func,
    post_trial_func = post_trial_func
    )

for rv in values(mcs_test.rvdict)
    @test rv.dist isa Mimi.SampleStore
end

for i = 1:num_scenarios
    @test results[i, :] == results[1, :]
end

rm(output_dir, recursive = true)

# Test in the case of sampling = LHS

generate_trials!(mcs_test, num_trials; sampling = LHS)
trial1 = copy(collect(values(mcs_test.rvdict))[1].dist.values)

for rv in values(mcs_test.rvdict)
    @test rv.dist isa Mimi.SampleStore
end

generate_trials!(mcs_test, num_trials; sampling = LHS)
trial2 = copy(collect(values(mcs_test.rvdict))[1].dist.values)

@test trial1!=trial2