text:
module HMDB

using LightXML
using Printf
import NMR
import Dates
import LinearAlgebra

global HMDB_dir
global hmdb_root

mutable struct HMDBpeaks
    spectID::String
    accession::String
    name::String
    nucleus::String
    frequency::String
    solvent::String
    pH::String
    pks::Array{Float64,1}
    ints::Array{Float64,1}
end

function HMDBpeaks(pks,ints; spectID="", accession="", name="", nucleus="",
    frequency="",solvent="",pH="")

    return HMDBpeaks(spectID,accession,name,nucleus,frequency,
      solvent,pH,pks,ints)
end

function HMDBpeaks(spectrum::XMLElement; accession="", name="", normalisation=1)
        peaks=find_element(spectrum,"nmr-one-d-peaks");
        snucleus=content(find_element(spectrum,"nucleus"));
        ssolvent=content(find_element(spectrum,"solvent"));
        sfreq=content(find_element(spectrum,"frequency"));
        spH=content(find_element(spectrum,"sample-ph"));
        sid=content(find_element(spectrum,"id"));
        pks=[parse(Float64,content(find_element(x,"chemical-shift"))) for x in child_elements(peaks)];
        ints=[parse(Float64,content(find_element(x,"intensity"))) for x in child_elements(peaks)];
        ints .*= normalisation/sum(ints)
	return HMDBpeaks(pks,ints,name=name,accession=accession,nucleus=snucleus,
                spectID=sid,frequency=sfreq,solvent=ssolvent,pH=spH)
end

global refPeaks

function __init__()
	global hmdb_root=root(parse_file("$(@__DIR__)/HMDB_subset_with_NMR.xml"))
	global refPeaks=Dict{String,HMDBpeaks}()
	name=""
	for m in hmdb_root["metabolite"]
		name = content(find_element(m,"name"))
		acc = content(find_element(m,"accession"))
		refspect = find_element(m,"nmr-ref-1D1H")
		spectProtons = parse(Int64, attribute(refspect,"spectProtons"))
		peaks = HMDBpeaks(find_element(refspect,"nmr-one-d"),accession=acc,normalisation=spectProtons,name=name)
		push!(refPeaks, name => peaks)
	end

	println("HMDB initialised from $(@__DIR__)");
end

"""
	(score,iscore,alpha,std) = matchPeaks(p::NMR.PeakStruct,ref::HMDBpeaks;tol=0.001)

match the peaks in `p` to the reference `ref`. A tuple of four values is returned:
- `score`: sum of the intensities of the peaks found, normalised by the height of the largest reference peak.  Missing peaks are counted as negative
- `iscore`: score normalised by the total intensity of the reference peaks (range -1...1)
- `alpha`: fitted normalised concentration
- `std`: mean error of fitted normalised concentration
"""
function matchPeaks(p::NMR.PeakStruct,ref::HMDBpeaks;tol=0.001)
    d=ref.pks
    i=ref.ints
    length(d) == length(i) || throw("lengths if peak position and intensity lists must match")

    npeaks=length(p.positions)
    score=0.0;

    q2=0.0; # cumulative normalisation
    p2=0.0;
    Pk=0.0; # cumulative projection

    for k=1:length(d)
        closest=p.positions[1];
        cInt=p.intensities[1];

        r=searchsorted(p.positions,d[k])
        if first(r) > npeaks
            closest=p.positions[npeaks]
            cInt=p.intensities[npeaks];

        elseif first(r)==1 && isempty(r)
            closest=p.positions[1]
            cInt=p.intensities[1]

        elseif !isempty(r)
            closest=p.positions[first(r)]
            cInt=p.intensities[first(r)]
        else
            closest= abs(p.positions[first(r)]-d[k])<abs(p.positions[first(r)-1]-d[k]) ? p.positions[first(r)] : p.positions[first(r)-1]
            cInt= abs(p.positions[first(r)]-d[k])<abs(p.positions[first(r)-1]-d[k]) ? p.intensities[first(r)] : p.intensities[first(r)-1]

        end

        if abs(closest - d[k]) < tol
          score += i[k]
            q2 += i[k]*i[k]
            Pk += i[k]*cInt
            p2 += cInt*cInt
        else
            score -= i[k]
        end
    end
    return (score,score/sum(i),Pk/q2, (p2-Pk*Pk/q2)/q2/(length(d)-1),ref.name)
end

function matchPeaks(p::NMR.PeakStruct,name::String; opts...)
	return matchPeaks(p, refPeaks[name]; opts...)
end

"""
	xml=matchReport(p::NMR.PeakStruct,tol=0.01,AutoShift=true,iscoreCutoff=0,Id="",verbose=true)

match peaks in `p` against all metabolites in `HMDB.refPeaks`, and create a table of the results.
If `AutoShift` is set to `true` (default), the `p` is automatically shfited for the greatest overlap
with D-Glucose. `tol` is the matching tolerance for peak identification. Metabolites
with iscores below `iscoreCutoff` are not reported. Metabolites are listed in decreasing
order of concentration. `matchReport()` returns an XML object containing the listed data. `Id` is incorporated
into the XML object as an attribute, to facilitate later identification of the match report.
"""
function matchReport(pinput::NMR.PeakStruct;tol=0.01,AutoShift=true,iscoreCutoff=0,verbose=true,Id="")
	xml=new_element("HMDBMatchReport");
	set_attribute(xml,"Id",Id);
	set_attribute(xml,"Date","$(Dates.now())")
	if AutoShift
		δ=-0.02:0.000125:0.02
		m=[HMDB.matchPeaks(NMR.shift(pinput,x),HMDB.refPeaks["D-Glucose"],tol=0.01)[2] for x in δ]
		p=NMR.shift(pinput,δ[argmax(m)])
	else
		p=pinput
	end

	matches=[HMDB.matchPeaks(p,ref,tol=tol) for (names,ref) in HMDB.refPeaks];
	filter!(x->x[2]>0,matches);         # remove all metabolites for which iscore<0
	sort!(matches,by=x->x[3],rev=true); # sort by descending concentration

	if verbose
		@printf("Metabolite Matches\n")
		@printf("==================\n")
		@printf("%-30s % 10s % 10s % 10s\n\n","Name","Score","iScore","alpha");
	end

	for (score,iscore,alpha,delta,name) in matches
   	 	verbose && @printf("%-30s % 10.2f% 10.2f% 14.3f ±%10.3f\n",name,score,iscore,alpha,sqrt(abs(delta)))
		c=new_child(xml,"Match")
		set_attributes(c,Name=name,score="$(score)", iscore="$(iscore)", alpha="$(alpha)",stderr="$(sqrt(abs(delta)))")
	end
	return xml
end

"""
    s::Data1D = refSpectrum(name::String,range;lw=0.01)

compute a reference spectrum over `range` using the peak positions in
reference spectrum `name`. A Lorentzian line with half width `lw` is used. The
reference spectrum is returned as a real `Data1D` object, unless the option `complex=true` is given.
"""
function refSpectrum(name::String,range::Array{Float64,1};lw=0.005,excl=x->false, complex=false)
	p=refPeaks[name]
  ints=[]
  pks=[]
  # filter out excluded peaks
  npeaks=0;
  for k=1:length(p.pks)
    if(!excl(p.pks[k]))
        push!(pks,p.pks[k])
        push!(ints,p.ints[k])
        npeaks+=1
      end
  end
    if (complex)
        vals = sum(k -> ints[k]*[NMR.clorentzian(pks[k],1.0/lw^2,x) for x in range], 1:npeaks)
        d=NMR.Data1D{Complex{Float64},Float64}(vals,minimum(range),maximum(range))
   else
        vals = sum(k -> ints[k]*[NMR.lorentzian(pks[k],1.0/lw^2,x) for x in range], 1:npeaks)
        d=NMR.Data1D{Float64,Float64}(vals,minimum(range),maximum(range))
   end
	return d
end

function refSpectrum(name::String,d::NMR.Data1D;opts...)
	return refSpectrum(name,collect(NMR.ind(d));opts...)
end

"""
	decompositionMatrix(names,range)

compute a decomposition matrix B such that B s  is the vector of concentrations
of the metabolites listed in `names`. The pseudo-inverse is used for this purpose.
"""
function decompositionMatrix(mlist,range;opts...)
	B=zeros(length(range),length(mlist))
	k=1
	for m in mlist
		r=HMDB.refSpectrum(m,range;opts...)
  	 	B[:,k]=NMR.val(r)
		k=k+1
	end
	return LinearAlgebra.pinv(B)
end

end

max_stars_repo_name: marcel-utz/NMR

max_issues_repo_name: marcel-utz/NMR


text:
using Plots, Measures; pyplot()

pDiscrete = [0.25, 0.25, 0.5]
xGridD = 0:2

pContinuous(x) = 3/4*(1 - x^2)
xGridC = -1:0.01:1

pContinuous2(x) = x < 0 ? x+1 : 1-x

p1 = plot(xGridD, line=:stem, pDiscrete, marker=:circle, c=:blue, ms=6, msw=0)
p2 = plot(xGridC, pContinuous.(xGridC), c=:blue)
p3 = plot(xGridC, pContinuous2.(xGridC), c=:blue)

plot(p1, p2, p3, layout=(1,3), legend=false, ylims=(0,1.1), xlabel="x",
	ylabel=["Probability" "Density" "Density"], size=(1200, 400), margin=5mm)

ext: jl

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
const SSH_MSG_DISCONNECT       = 1
const SSH_MSG_IGNORE           = 2
const SSH_MSG_UNIMPLEMENTED    = 3
const SSH_MSG_DEBUG            = 4
const SSH_MSG_SERVICE_REQUEST  = 5
const SSH_MSG_SERVICE_ACCEPT   = 6
const SSH_MSG_KEXINIT          = 20
const SSH_MSG_NEWKEYS          = 21
const SSH_MSG_KEXDH_INIT       = 30
const SSH_MSG_KEXDH_REPLY      = 31
const SSH_MSG_USERAUTH_REQUEST = 50
const SSH_MSG_USERAUTH_FAILURE = 51
const SSH_MSG_USERAUTH_SUCCESS = 52
const SSH_MSG_USERAUTH_BANNER  = 53

# Messages for Authentication method publickey
const SSH_MSG_USERAUTH_PK_OK   = 60

# Messages for Authentication method keyboard-interactive
const SSH_MSG_USERAUTH_INFO_REQUEST     =  60
const SSH_MSG_USERAUTH_INFO_RESPONSE    =  61

# Messages for ssh-connection
const SSH_MSG_GLOBAL_REQUEST            =  80
const SSH_MSG_REQUEST_SUCCESS           =  81
const SSH_MSG_REQUEST_FAILURE           =  82
const SSH_MSG_CHANNEL_OPEN              =  90
const SSH_MSG_CHANNEL_OPEN_CONFIRMATION =  91
const SSH_MSG_CHANNEL_OPEN_FAILURE      =  92
const SSH_MSG_CHANNEL_WINDOW_ADJUST     =  93
const SSH_MSG_CHANNEL_DATA              =  94
const SSH_MSG_CHANNEL_EXTENDED_DATA     =  95
const SSH_MSG_CHANNEL_EOF               =  96
const SSH_MSG_CHANNEL_CLOSE             =  97
const SSH_MSG_CHANNEL_REQUEST           =  98
const SSH_MSG_CHANNEL_SUCCESS           =  99
const SSH_MSG_CHANNEL_FAILURE           = 100

const SSH_OPEN_ADMINISTRATIVELY_PROHIBITED = 1
const SSH_OPEN_CONNECT_FAILED              = 2
const SSH_OPEN_UNKNOWN_CHANNEL_TYPE        = 3
const SSH_OPEN_RESOURCE_SHORTAGE           = 4

const SSH_DISCONNECT_HOST_NOT_ALLOWED_TO_CONNECT    =  1
const SSH_DISCONNECT_PROTOCOL_ERROR                 =  2
const SSH_DISCONNECT_KEY_EXCHANGE_FAILED            =  3
const SSH_DISCONNECT_RESERVED                       =  4
const SSH_DISCONNECT_MAC_ERROR                      =  5
const SSH_DISCONNECT_COMPRESSION_ERROR              =  6
const SSH_DISCONNECT_SERVICE_NOT_AVAILABLE          =  7
const SSH_DISCONNECT_PROTOCOL_VERSION_NOT_SUPPORTED =  8
const SSH_DISCONNECT_HOST_KEY_NOT_VERIFIABLE        =  9
const SSH_DISCONNECT_CONNECTION_LOST                = 10
const SSH_DISCONNECT_BY_APPLICATION                 = 11
const SSH_DISCONNECT_TOO_MANY_CONNECTIONS           = 12
const SSH_DISCONNECT_AUTH_CANCELLED_BY_USER         = 13
const SSH_DISCONNECT_NO_MORE_AUTH_METHODS_AVAILABLE = 14
const SSH_DISCONNECT_ILLEGAL_USER_NAME              = 15


text:
module BenchNondeterministic

using BenchmarkTools
using Folds
using FoldsThreads

partially_parallelizable(seq) = (gcd(y, 42) for x in seq for y in 1:10000x)

demo(ex) = Folds.sum(partially_parallelizable(Iterators.Stateful(1:100)), ex)

const SUITE = BenchmarkGroup()
SUITE["seq"] = @benchmarkable(demo(SequentialEx()))
SUITE["nondet"] = @benchmarkable(demo(NondeterministicEx()))

end  # module
BenchNondeterministic.SUITE


text:
module PyUtils
    using Reexport
    @reexport using PyCall

    export jl_to_np_type, np_to_jl_type
    export PyArray_ReadOnly

    const _pyslice = PyNULL()
    const _pycolon = PyNULL()

    # Upstream fixes and new conversions
    include("PyCall.jl")

    include("scalar_type_conversions.jl")

    function __init__()
        copy!(_pyslice, py"slice")
        copy!(_pycolon, _pyslice(nothing, nothing, nothing))
    end
end


text:
abstract type AbstractWalker end

# Make walker structs
struct Walker1D{T} <: AbstractWalker
    x::T
end

struct Walker2D{T} <: AbstractWalker
    x::T
    y::T
end

struct Step2D{T} <: AbstractWalker
    x::T
    y::T
end

struct Step1D{T} <: AbstractWalker
    x::T
end

struct walker_end_point_meta <: AbstractWalker
    time::Float64
    rate::Float64
    xₜ::Int64
    yₜ::Int64
end

struct WalkersData <:AbstractWalker
    nwalks::Int
    timesteps::Int
    rate::Float64
    width::Int
    length::Int
    x::Int
    y::Int
end


text:
##################################################################################
# AFFINE TRANSITION
##################################################################################

"""
Convert a tuple of coordinates to a "vector" format.

# Arguments
- `tuple`::Tuple{Float, Float}`: point coordinate (y, x).
"""
function tupleToVector(tuple)
    return [elem for elem in tuple]
end

"""
Create a transition matrix for a series of operations.
Each parameter in input list defines a type (scaling or rotation) of operation and its ratio/angle.

# Arguments
- `parameters`::Array{Tuple{Char, Float/Tuple{Float, Float}}, 1}`: point coordinate (y, x).
"""
function mapMatrix(parameters)
    matrix = [1 0; 0 1]
    for trans in parameters
        if trans[1] == 's'
            matrix = [trans[2][1] 0; 0 trans[2][2]] * matrix
        elseif trans[1] == 'r'
            angleRad = - trans[2] / 180 * pi
            matrix = [cos(angleRad) -sin(angleRad); sin(angleRad) cos(angleRad)] * matrix
        end
    end
    return matrix
end

"""
Find coordinates of a point after transition by a given matrix.

# Arguments
- `y`::Float` y coordinate of a point.
- `x`::Float` x coordinate of a point.
- `originCnvrtd`::Array{Float, 1}` origin point as a "vector".
- `mapMtrx`::Array{Float, 2}` transition matrix.
"""
function mappedCoords(y, x, originCnvrtd, mapMtrx)
    return reverse(mapMtrx * ([x; y] .- reverse(originCnvrtd)) .+ reverse(originCnvrtd))
end

"""
Find ranges of selection area after a transition by a matrix.
Available both transformation and its inversion evaluation.

# Arguments
- `colorMatrix`::Array{Float, 2}` full-size matrix of one color.
- `ld::Tuple{Integer, Integer}`: lower-left selection area corner coordinates as (y, x).
- `ru::Tuple{Integer, Integer}`: upper-right selection area corner coordinates as (y, x).
- `origin::Tuple{Number, Number}`: transition origin point as (y, x).
- `mapMtrx`::Array{Float, 2}` transition matrix.
- `translVector`::Array{Float, 1}` translation vector to move floating area by.
- `rever`::Bool` information about straight or reversed checking.
"""
function tRange(colorMatrix, ld, ru, origin, mapMtrx, translVector, rever=false)
    sizeOfM = size(colorMatrix)
    crnrs = tupleToVector.([ld, ru, (ld[1], ru[2]), (ru[1], ld[2])])
    crnrsMapped = [mappedCoords(corner[1] - translVector[1] * (rever), corner[2] - translVector[2] * (rever), origin, mapMtrx) .+ translVector .* (!rever) for corner in crnrs]
    choiceY = [corner[1] for corner in crnrsMapped]
    choiceX = [corner[2] for corner in crnrsMapped]

    boundsPreY = Int.([floor(minimum(choiceY)), ceil(maximum(choiceY))])
    boundsPreX = Int.([floor(minimum(choiceX)), ceil(maximum(choiceX))])
    if rever == false
        ldNew = (minimum([sizeOfM[1], boundsPreY[2]]), maximum([1, boundsPreX[1]]))
        ruNew = (maximum([1, boundsPreY[1]]), minimum([sizeOfM[2], boundsPreX[2]]))

        return [[ruNew[1], ldNew[1]], [ldNew[2], ruNew[2]]]
    else
        ldNew = (boundsPreY[2], boundsPreX[1])
        ruNew = (boundsPreY[1], boundsPreX[2])
        return [[ruNew[1], ldNew[1]], [ldNew[2], ruNew[2]]]
    end
end

"""
Find the extra frame dimensions.
They are used to avoid trying to evaluate pixels from out of a initial selection range.

# Arguments
- `tRepRangeBnds1`::Array{Integer, 2}` vertical boundaries after inverted ranges evaluation.
- `tRepRangeBnds2`::Array{Integer, 2}` horizontal boundaries after inverted ranges evaluation.
- `ld::Tuple{Integer, Integer}`: lower-left selection area corner coordinates as (y, x).
- `ru::Tuple{Integer, Integer}`: upper-right selection area corner coordinates as (y, x).
"""
function mtrxExtensionRngs(tRepRangeBnds1, tRepRangeBnds2, ld, ru)
    extRanges = [[max(ru[1] - tRepRangeBnds1[1] - 1, 0), max(tRepRangeBnds1[2] - ld[1] - 1, 0) + 1],
        [max(ld[2] - tRepRangeBnds2[1] - 1, 0), max(tRepRangeBnds2[2] - ru[2] - 1, 0) + 1]]
    return extRanges
end

"""
Prepare matrices of selections extented by a proper frame.

# Arguments
- `extRngs`::Array{Array{Integer, 1}, 1}` vertical boundaries after inverted ranges evaluation.
- `colorMatrix`::Array{Float, 2}` selection from a matrix of one color.
- `mode::Integer` color mode is 3; mask (artificial-alpha) mode is 1.
"""
function extendedInitialMatrix(extRngs, colorMatrix, mode)
    noneValue = 1
    if mode == 1
        edgeValue = 1
        sizeYOfM = size(colorMatrix, 1)
        horizontal = hcat(fill(noneValue, (sizeYOfM, extRngs[2][1])), fill(edgeValue, (sizeYOfM, 1)),
        colorMatrix, fill(edgeValue, (sizeYOfM, 1)), fill(noneValue, (sizeYOfM, extRngs[2][2])))
        sizeXOfH = size(horizontal, 2)
        return vcat(fill(noneValue, (extRngs[1][1], sizeXOfH)), fill(edgeValue, (1, sizeXOfH)),
            horizontal, fill(edgeValue, (1, sizeXOfH)), fill(noneValue, (extRngs[1][2], sizeXOfH)))

    else
        sizeYOfM = size(colorMatrix, 1)
        horizontal = hcat(fill(noneValue, (sizeYOfM, extRngs[2][1])), colorMatrix[1:end, 1],
        colorMatrix, colorMatrix[1:end, end], fill(noneValue, (sizeYOfM, extRngs[2][2])))
        sizeXOfH = size(horizontal, 2)
        return vcat(fill(noneValue, (extRngs[1][1], sizeXOfH)), horizontal[1, 1:end]',
            horizontal, horizontal[end, 1:end]', fill(noneValue, (extRngs[1][2], sizeXOfH)))
    end
end

"""
Calculate (by bilinear interpolation) the value for pixel which coordinates of origin are given.
Theese coordinates are not integer, so the calculation is inevitable.

# Arguments
- `y`::Float` y coordinate of a point.
- `x`::Float` x coordinate of a point.
- `colorMatrix`::Array{Float, 2}` (extended) selection from a matrix of one color.
"""
function evaluate(y, x, colorMatrix)
    xFloor = Int(floor(x))
    yFloor = Int(floor(y))

    difs = [[1 - y + yFloor, -yFloor + y], [1 - x + xFloor, -xFloor + x]]
    return sum([difs[1][1] * difs[2][1] * colorMatrix[yFloor, xFloor],
                difs[1][1] * difs[2][2] * colorMatrix[yFloor, xFloor + 1],
                difs[1][2] * difs[2][1] * colorMatrix[yFloor + 1, xFloor],
                difs[1][2] * difs[2][2] * colorMatrix[yFloor + 1, xFloor + 1]])
end

"""
Get the selection and its position after a series of given transformations on RGB image.

# Arguments
- `initValues::Array{Array{Float, 2}, 1}` series of input color or mask matrices in.
- `ld::Tuple{Integer, Integer}`: lower-left selection area corner coordinates as (y, x).
- `ru::Tuple{Integer, Integer}`: upper-right selection area corner coordinates as (y, x).
- `origin::Tuple{Number, Number}`: transition origin point as (y, x).
- `parameters::Array{Tuple{Char, Float/Tuple{Float, Float}}, 1}`: point coordinate (y, x).
- `translVectorTpl`::Tuple{Float, Float}` raw translation vector to move floating area by.
"""
function transform(initValues, ld, ru, origin, parameters, translVectorTpl)
    translVector = [translVectorTpl[1]; translVectorTpl[2]]
    final =  Array{Array{Float64,2},1}(undef, length(initValues))

    mapMtrx = mapMatrix(parameters)
    invMapMtrx = zeros(2, 2)
    try
        invMapMtrx = inv(mapMtrx)
    catch
        throw(ErrorException("Cannot scale by 0 value."))
    end

    originCnvrtd = tupleToVector(origin)

    acqrRanges = tRange(initValues[1], ld, ru, origin, mapMtrx, translVector)
    tRange1 = acqrRanges[1][1]:acqrRanges[1][2]
    tRange2 = acqrRanges[2][1]:acqrRanges[2][2]
    repLd = (acqrRanges[1][2], acqrRanges[2][1])
    repRu = (acqrRanges[1][1], acqrRanges[2][2])
    acqrRepRanges = tRange(initValues[1][ru[1]:ld[1], ld[2]:ru[2]], repLd, repRu, origin, invMapMtrx, translVector, true)
    extRngs = mtrxExtensionRngs([acqrRepRanges[1][1], acqrRepRanges[1][2]],
             [acqrRepRanges[2][1], acqrRepRanges[2][2]], ld, ru)

    if length(tRange1) == 0 || length(tRange2) == 0
        @simd  for col in 1:length(initValues)
            final[col] = ones(size(initValues[col]))
        end
        return (final, ((0, 0), (0, 0)))

    else
        @simd  for col in 1:length(initValues)
            extndMatrix = extendedInitialMatrix(extRngs, initValues[col][ru[1]:ld[1], ld[2]:ru[2]], length(initValues))
            sectTransformed = [evaluate((mappedCoords(i - translVector[1], j - translVector[2], originCnvrtd, invMapMtrx)
            .+ [extRngs[1][1] + 1, extRngs[2][1] + 1] .- [ru[1] - 1, ld[2] - 1])..., extndMatrix) for i in tRange1, j in tRange2]
            final[col] = ones(size(initValues[col]))
            final[col][tRange1, tRange2] = sectTransformed
        end
    end
    return (final, ((tRange1[end], tRange2[1]), (tRange1[1], tRange2[end])))
end

##################################################################################
# GENERAL TRANSFORMATION FUNCTIONS
##################################################################################

"""
Give a deafult origin point - at the center of selection.

# Arguments
- `ld::Tuple{Integer, Integer}`: lower-left scaling area corner coordinates as (y, x).
- `ru::Tuple{Integer, Integer}`: upper-right scaling area corner coordinates as (y, x).
"""
function defaultOrigin(ld, ru)
    return (ru[1]+(ld[1]-ru[1])/2, ld[2]+(ru[2]-ld[2])/2)
end

"""
Check the validity of corner values given on input. Return true if correct.

# Arguments
- `initValues::Array{Array{Float64,2},1}`: r, g, b matrices list.
- `ld::Tuple{Integer, Integer}`: lower-left scaling area corner coordinates as (y, x).
- `ru::Tuple{Integer, Integer}`: upper-right scaling area corner coordinates as (y, x).
"""
function validity(initValues, ld, ru)
    vals = [ld..., ru...]
    if ! all((typeof.(vals) ) .<: Number)
        return false
    elseif any(vals .< 1)
    return false
  elseif ld[1] > size(initValues[1], 1) || ru[2] > size(initValues[1], 2)
    return false
  elseif ld[1] < ru[1] || ld[2] > ru[2]
    return false
  else
    return true
    end
end

"""
Paste transformed selection into the image.

# Arguments
- `colorMatrices::Array{Array{Float, 2}, 1}`: previous color matrices.
- `transformedMatrices::Array{Array{Float, 1}, 1}`: image selection after a series of transformation.
- `transformedMask::Array{Array{Float, 1}, 1}`: mask selection after a series of transformation.
- `ld::Tuple{Integer, Integer}`: lower-left scaling area corner coordinates as (y, x).
- `ru::Tuple{Integer, Integer}`: upper-right scaling area corner coordinates as (y, x).
"""
function merge(colorMatrices, transformedMatrices, transformedMask, ld, ru)
    outImage = Array{Array{Float64,2},1}(undef, 3)
    @simd for col in 1:3
        outImage[col] = colorMatrices[col]
        outImage[col][ru[1]:ld[1], ld[2]:ru[2]] .= 1
        outImage[col] =  transformedMatrices[col] .* (1 .- transformedMask[1]) + transformedMask[1] .* outImage[col]
    end
    return outImage
end

"""
Conduct transformations.

# Arguments
- `initValues::Array{Array{Float, 2}, 1}`: previous color matrices.
- `ld::Tuple{Integer, Integer}`: lower-left scaling area corner coordinates as (y, x).
- `ru::Tuple{Integer, Integer}`: upper-right scaling area corner coordinates as (y, x).
- `transfList::Array{Tuple{Tuple{Float, Float}, /parameters/, Tuple{Float, Float}}, 1}`:
array of operations on image, each containing origin, parameters(described previously), translation vector.

# Example
selectionTransform(initValues, (400, 300), (200, 700),
 [ ((300, 500), [('r', 100), ('s', (2, 1))], (20.4, -40)) ])

Does the transformations on initValues matrices, where selecion is defined by LD = (400, 300), RU = (200, 700)
and the origin point is (300, 500). Image is translated by (20.4, -40), rotated by 100 degrees and scaled with ratios 2, 1.
"""
function selectionTransform(initValues, ld, ru, transfList)
    if ! validity(initValues, ld, ru)
        throw(ErrorException("Selection range error. Check LD and RU coordinates."))
    end
    mask = [zeros(size(initValues[1]))]
    floating = initValues
    floatingLd = ld
    floatingRu = ru
    for operation in transfList
        floating = transform(floating, floatingLd, floatingRu, operation[1], operation[2], operation[3])[1]
        tempData = transform(mask, floatingLd, floatingRu, operation[1], operation[2], operation[3])
        mask = tempData[1]
        floatingLd = tempData[2][1]
        floatingRu = tempData[2][2]
        if floatingLd == (0, 0) || floatingRu == (0, 0)
            throw(ErrorException("Cannot transform by given parameters. Range error."))
        end
    end
    return merge(initValues, floating, mask, ld, ru)
end

##################################################################################
# SYMMETRY
##################################################################################
function mirror(RGB, ax)
    if ax == :x
        dim = 1
    elseif ax == :y
        dim = 2
    else
        return
    end
    newRGB = reverse.(RGB, dims=dim)
    return newRGB
end

max_stars_repo_name: o-Mateo-o/ImageEditorApp

max_issues_repo_name: o-Mateo-o/ImageEditorApp

max_forks_repo_name: o-Mateo-o/ImageEditorApp


text:
# TODO: Use LinuxPerf.jl after https://github.com/JuliaPerf/LinuxPerf.jl/pull/22

struct PerfConfig
    detailed::Int
    event::Union{Nothing,String}
    per_thread::Bool
    all_cpus::Bool
    cpu::Union{Nothing,String}
    _perfopts::Cmd  # (for easily trying new options) TODO: remove this
end

function as_perf_config(;
    detailed::Union{Integer,Nothing} = nothing,
    event::Union{AbstractVector{<:AbstractString},Nothing} = nothing,
    per_thread::Bool = false,
    cpu::Union{AbstractString,Nothing} = nothing,
    all_cpus::Bool = cpu !== nothing,
    _perfopts::Cmd = ``,
)
    if detailed === nothing
        detailed = 0
    elseif detailed < 0
        error("`detailed` must be 1, 2, or 3")
    end
    if event !== nothing
        event = join(event, ",")
    end
    return PerfConfig(detailed, event, per_thread, all_cpus, cpu, _perfopts)
end

const CacheDict = Dict{String,Union{Float64,Int}}
const CacheRef = Ref{Union{Nothing,CacheDict}}

struct TrialPerf
    output::String
    cache::CacheRef
end

TrialPerf(output) = TrialPerf(output, CacheRef(nothing))
JSON.lower(t::TrialPerf) = (; output = getfield(t, :output))

struct GroupPerf
    group::Dict{Any,Union{TrialPerf,GroupPerf}}
end

const GenericPerf = Union{TrialPerf,GroupPerf}

GroupPerf() = GroupPerf(Dict{Any,GenericPerf}())

Base.keys(g::GroupPerf) = keys(g.group)
Base.get(g::GroupPerf, k, default) = get(g.group, k, default)
Base.getindex(g::GroupPerf, k) = g.group[k]
Base.setindex!(g::GroupPerf, v::GenericPerf, k) = g.group[k] = v

const FIELD_SEPARATOR = '\\'  # as recommended in perf-stat(1)

function Base.Cmd(cfg::PerfConfig)
    cmd = `perf stat`
    if cfg.detailed > 0
        detailed = "-" * "d"^cfg.detailed
        cmd = `$cmd $detailed`
    end
    if cfg.event !== nothing
        cmd = `$cmd --event=$(cfg.event)`
    end
    if cfg.per_thread
        cmd = `$cmd --per-thread`
    end
    cmd = `$cmd $(cfg._perfopts)`
    cmd = `$cmd --field-separator=$FIELD_SEPARATOR`
    if cfg.all_cpus
        cmd = `$cmd --all-cpus`
        if cfg.cpu !== nothing
            cmd = `$cmd --cpu=$(cfg.cpu)`
        end
    else
        cmd = `$cmd --pid=$(getpid())`
    end
    return cmd
end

function withperf(f, cfg::PerfConfig)
    cmd = Cmd(cfg)
    output = Ref{String}()
    outpipe = Pipe()
    proc = run(pipeline(cmd, stdout = stderr, stderr = outpipe); wait = false)
    local y
    @sync begin
        try
            close(outpipe.in)
            @async output[] = read(outpipe, String)
            y = f()
        catch
            close(outpipe)
            rethrow()
        finally
            @debug "Stopping perf:" cmd
            flush(stderr)
            kill(proc, Base.SIGINT)
            wait(proc)
        end
    end
    return (y, TrialPerf(output[]))
end

function getcache(trial::TrialPerf)
    ref = getfield(trial, :cache)
    cache = ref[]
    if cache === nothing
        ref[] = cache = parsestat(getfield(trial, :output))
    end
    return cache
end

function parsestat(output)
    cache = CacheDict()
    for ln in eachline(IOBuffer(output))
        parts = split(ln, FIELD_SEPARATOR)
        isempty(parts[1]) && continue
        k = parts[3]
        parts[1] == "<not supported>" && continue
        v = tryparse(Int, parts[1])
        if v === nothing
            v = tryparse(Float64, parts[1])
            if v === nothing
                error("failed parsing row: ", ln)
            end
        end
        cache[k] = v
    end
    return cache
end

Base.keys(trial::TrialPerf) = keys(getcache(trial))
Base.get(trial::TrialPerf, key, default) = get(getcache(trial), key, default)
Base.getindex(trial::TrialPerf, event::AbstractString) = getcache(trial)[event]

function Base.propertynames(::TrialPerf, private::Bool = false)
    pubnames = (ALL_NAMES..., :ratio, :percent)
    if private
        return (pubnames..., fieldnames(TrialPerf))
    else
        return pubnames
    end
end

function Base.getproperty(trial::TrialPerf, name::Symbol)
    if name in ALL_NAMES
        return getvalue(trial, Val{name}())
    elseif name === :ratio
        return RatioAccessor(trial)
    elseif name === :percent
        return PercentAccessor(trial)
    else
        return getfield(trial, name)
    end
end

fulltable(trial::TrialPerf) = fulltable(trial.output)
function fulltable(output::AbstractString)
    table = (
        value = Union{Float64,Int64,Nothing}[],
        unit = Union{String,Nothing}[],
        event = String[],
        runtime = Union{Int64,Nothing}[],
        percentage = Union{Float64,Nothing}[],
        optional_metric = Union{Float64,Nothing}[],
        optional_unit = Union{String,Nothing}[],
    )
    for ln in eachline(IOBuffer(output))
        parts = split(ln, FIELD_SEPARATOR)
        push!(
            table.value,
            something(tryparse(Int, parts[1]), tryparse(Float64, parts[1]), Some(nothing)),
        )
        push!(table.unit, isempty(parts[2]) ? nothing : parts[2])
        push!(table.event, parts[3])
        push!(table.runtime, tryparse(Int64, parts[4]))
        push!(table.percentage, tryparse(Float64, parts[5]))
        push!(table.optional_metric, tryparse(Float64, parts[6]))
        push!(table.optional_unit, parts[7])
    end
    return table
end


text:
#---------------------------------------------------------------------#
#This function computes the global Mass and Differentiation Matrices.
#Written by F.X. Giraldo on April 19 2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function create_rhs(q,u,Dhat,Fhat,intma,DFloat)

    #Get lengths of arrays
    (Np,nelem)=size(intma)
    Npoin=size(q)

    #store flux term
    f=q*u

    #build numerical flux: fstar
    fstar=zeros(DFloat,Npoin)
    for e=1:nelem
        L=intma[Np,e]
        ep1=e+1
        if (e == nelem)
            ep1=1
        end
        R=intma[1,ep1]
        flux=0.5*( f[L] + f[R] - u*(q[R] - q[L]) )
        fstar[L]=flux
        fstar[R]=flux
    end

    #Build RHS vector
    R=Dhat*f - Fhat*fstar

    #Return arrays
    return(R)
end


text:
using Documenter, NetCDF

makedocs(
    modules = [NetCDF],
    clean   = false,
    format   = :html,
    sitename = "NetCDF.jl",
    authors = "Fabian Gans and contributors",
    pages = [
      "Home" => "index.md",
      "Manual" => [
        "quickstart.md",
        "highlevel.md",
        "intermediate.md",
        "strings.md"
      ]
    ]
)

deploydocs(
    #deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo   = "github.com/JuliaGeo/NetCDF.jl.git",
    julia  = "1.0",
    deps   = nothing,
    make   = nothing,
    target = "build"
)


text:
"""
        function get_rescaled_contact_matrices(filename::AbstractString)

Get three age-group contact matrices: contacts at home, work and social contacts, and, contacts
at school. The `JLD2` file at `filename` must contain the underlying contact matrices.
"""
function get_rescaled_contact_matrices(filename::AbstractString)
        M_Data = FileIO.load(filename)
        M_Kenya = M_Data["M_Kenya"]
        M_Kenya_ho = M_Data["M_Kenya_ho"]
        M_Kenya_other = M_Data["M_Kenya_other"]
        M_Kenya_work = M_Data["M_Kenya_work"]
        M_Kenya_school = M_Data["M_Kenya_school"]
        Ronescaler = 1/Real(eigvals(M_Kenya)[end])
        M_Kenya_ho .= M_Kenya_ho.*Ronescaler
        M_Kenya_other .= M_Kenya_other.*Ronescaler
        M_Kenya_work .= M_Kenya_work.*Ronescaler
        M_Kenya_school .= M_Kenya_school.*Ronescaler
        return M_Kenya_ho,M_Kenya_other,M_Kenya_school,M_Kenya_work
end

"""
        function get_population_size_matrix(filename)

Get a `NamedArray` data structure for numbers in each age group and each county of Kenya, as per 2019 census.
"""
function get_population_size_matrix(filename)
        df_pop = DataFrame(CSV.File(filename,stringtype=String))
        #df_pop = DataFrame(CSV.File(filename))
        N = NamedArray([Float64(df_pop[i, a+1][1]) for a = 1:17, i = 1:length(df_pop.county)])
        setnames!(N, vcat([string((a - 1) * 5) * "-" * string(a * 5 - 1) for a = 1:16], ["80+"]), 1)
        setnames!(N, df_pop.county, 2)
        setdimnames!(N, ["Age", "County"])
        return N
end

"""
 Function to reduce the age catorgories from 17 to 6
 This function collapses age groups into [0-19], [20-49], [50-59], [60-69], [70-79], 80+
"""

function reduce_age_categories(M_age,N_kenya,county)
    N_county = N_kenya[:,county]
    P_b =  N_county/sum(N_county)
    p_b_prime = vcat(sum(P_b[1:4]),sum(P_b[5:10]),sum(P_b[11:12]),sum(P_b[13:14]),sum(P_b[15:16]),P_b[17])
    p_vec = [[1,2,3,4],[5,6,7,8,9,10],[11,12],[13,14],[15,16],[17]]
    M_age_prime = zeros(6,6)

    for a_prime in 1:6, b_prime in 1:6
        M_age_prime[a_prime,b_prime] = sum([M_age[a,b]*P_b[b]/p_b_prime[b_prime] for a in p_vec[a_prime] for b in p_vec[b_prime]])
    end
    return  M_age_prime
end


text:
# This file is a part of SimilaritySearch.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

export LogSatNeighborhood

struct LogSatNeighborhood <: NeighborhoodAlgorithm
    base::Float64
end

function LogSatNeighborhood()
    return LogSatNeighborhood(1.1)
end

function optimize_neighborhood!(algo::LogSatNeighborhood, index::SearchGraph{T}, dist::Function, perf, recall) where {T}
    # optimize_neighborhood!(algo.g, index, perf, recall)
end

function neighborhood(algo::LogSatNeighborhood, index::SearchGraph{T}, dist::Function, item::T) where {T}
    n = length(index.db)
    k = max(1, ceil(Int, log(algo.base, n)))
    knn = search(index, dist, item, KnnResult(k))
    N = Vector{Int32}(undef, 0)

    @inbounds for p in knn
        dqp = p.dist
        pobj = index.db[p.objID]
        near = KnnResult(1)
        push!(near, p.objID, p.dist)
        for nearID in N
            d = convert(Float32, dist(index.db[nearID], pobj))
            push!(near, nearID, d)
        end

        f = first(near)
        if f.objID == p.objID
            push!(N, p.objID)
        end
    end

    return knn, N
end


text:
using BrainForth: @bf, @run

@bf promptr = [["Enter your name: ", prompt], whileempty]
@bf greeting = ["Hello, ", print, print, "!", println]

@bf main = [promptr, dupv, greeting,
            "Here's your name in reverse: ", print,
            reverse, println]

@run [main]

max_stars_repo_name: JuliaTagBot/BrainForth

max_issues_repo_name: JuliaTagBot/BrainForth

max_forks_repo_name: JuliaTagBot/BrainForth


text:
function samplesystem(;name, τ=3.0)
    @variables t
    D = Differential(t)
    sts = @variables x(t) RHS(t)  # RHS is observed
    ps = @parameters τ=τ       # parameters
    ODESystem([ RHS  ~ (1 - x)/τ, D(x) ~ RHS ], t, sts, ps; name)
end

max_stars_repo_name: bgctw/MTKHelpers

max_issues_repo_name: bgctw/MTKHelpers

max_forks_repo_name: bgctw/MTKHelpers


text:
module PlotUtils

export get_custom_rdbu_scale, get_plot_resolution, map_bool_to_color

"""
	get_custom_rdbu_scale(0.75)

Returns a customized red blue scale adjusted for a very thin white colorband
in the middle of the interval.
"""
function get_custom_rdbu_scale(opacity::Union{Float32, Float64})
	@assert opacity >= 0 || opacity <= 1 "Opacity must be in [0, 1]"
	return [
		[0, "rgba(255, 0, 0, $opacity)"],
		[0.45, "rgb(255, 200, 100, $opacity)"],
		[0.5, "rgba(255, 255, 255, $opacity)"],
		[0.55, "rgb(100, 200, 255, $opacity)"],
		[1, "rgba(0, 0, 255, $opacity)"]
	]
end

"""
	get_plot_resolution("large")

Returns a tuple of a resolution with 16:9 aspect ratio.
"""
function get_plot_resolution(size::String)
	size = lowercase(size)
	@assert size in ["small", "medium", "large", "xlarge"] "Requested plot size '$size' is not supported."

	if size == "small"
		# 16:9 360p
		return 640, 360
	elseif size == "medium"
		# 16:9 480p
		return 854, 480
	elseif size == "large"
		# 16:9 720p
		return 1280, 720
	elseif size == "xlarge"
		# 16:9 1080p
		return 1920, 1080
	end
end

"""
	map_bool_to_color([false, true, true, false], "red", "blue")

Maps a boolean vector to the two specified colors.
"""
function map_bool_to_color(bool_vec::Union{BitVector, Vector{Bool}}, color_true::String, color_false::String)
	return map(el -> el == true ? color_true : color_false, bool_vec)
end

end


text:
@typeparams struct MutableBinaryHeapWithReduction
    v::{}
    r::{}
    k2i::{}
    i2k::{}
    order::{<:Ordering}
    reduce::{}
end

function MutableBinaryHeapWithReduction(
    v,
    r = Vector{eltype(v)}(undef, length(v)),
    i2k = nothing,
    k2i = nothing;
    reduce,
    lt = isless,
    by = identity,
    rev::Bool = false,
    order::Ordering = Forward
)
    if isnothing(i2k)
        @assert isnothing(k2i)
        n = length(v)
        i2k = collect(1:n)
        k2i = collect(1:n)
    elseif isnothing(k2i)
        k2i = Dict(k=>i for (i,k) in enumerate(i2k))
    end
    order = ord(lt, by, rev, order)
    heapify!(
        PermutationTracker(v,k2i,i2k),
        By(((k,vk),) -> vk, order)
    )
    reduce_subtrees!(v,r,reduce)
    return MutableBinaryHeapWithReduction(
        v, r, k2i, i2k,
        order, reduce
    )
end

function update!(h::MutableBinaryHeapWithReduction, k, vk)
    @unpack v,r,k2i,i2k,order,reduce = h
    n = length(v)
    if n == 0; return k; end

    i = k2i[k]
    tv = PermutationTracker(v, k2i, i2k)
    tvk = (k,vk)
    torder = By(((k,vk),) -> vk, order)

    if lt(order, vk, v[i])
        i = walk(
            piggyback(
                Bubbler!(tv, tvk, torder, Val(:up)),
                Reducer!(v,r,reduce)
            ),
            n, i, Val(:up)
        )
    else
        i = walk(
            Bubbler!(tv, tvk, torder, Val(:down)),
            n, i, Val(:down)
        )
    end
    tv[i] = tvk
    walk(Reducer!(v,r,reduce),n,i, Val(:up))
    return i
end


text:
using Flux
using StatsBase
using Statistics
using Distributions
using DistributionsAD
using SumProductTransform

struct SPTN
	m
end

Flux.@functor SPTN

function SPTN(;idim::Int=1, activation=identity, ncomp=4, nlayers::Int=2,
							init_seed=nothing, unitary=:butterfly, sharing=:dense, firstdense=false, kwargs...)
	# if seed is given, set it
	(init_seed != nothing) ? Random.seed!(init_seed) : nothing

	model = SumProductTransform.buildmixture(
				idim,
				ncomp,
				nlayers,
				activation;
				sharing = sharing,
				firstdense = firstdense,
				unitary = unitary)

	# reset seed
	(init_seed != nothing) ? Random.seed!() : nothing

	SPTN(model)
end

function Base.show(io::IO, sptn::SPTN)
	# to avoid the show explosion
	print(io, "SPTN(...)")
end

function StatsBase.predict(model::SPTN, x)
	-logpdf(model.m, x)
end

function StatsBase.fit!(model::SPTN, data::Tuple; max_train_time=82800,
						batchsize=64, patience=20, check_interval::Int=10, kwargs...)
	opt = ADAM()
	history = MVHistory()
	tr_model = deepcopy(model)
	ps = Flux.params(tr_model)
	_patience = patience

	# split data
	tr_x = data[1][1]
	val_x = data[2][1][:, data[2][2] .== 0]

	best_val_loss = Inf
	i = 1
	start_time = time()
	frmt = (v) -> round(v, digits=4)
	for batch in RandomBatches(tr_x, batchsize)
		# batch loss
		batch_loss = 0f0
		grad_time = @elapsed begin
			gs = gradient(() -> begin
				batch_loss = -mean(logpdf(tr_model.m, batch))
			end, ps)
			Flux.update!(opt, ps, gs)
		end

		push!(history, :training_loss, i, batch_loss)
		push!(history, :grad_time, i, grad_time)

		# validation/early stopping
		if (i%check_interval == 0)
			val_lkl_time = @elapsed val_loss = -mean(batchlogpdf(tr_model.m, val_x, batchsize))
			@info "$i - loss: $(frmt(batch_loss)) (batch) | $(frmt(val_loss)) (validation) || $(frmt(grad_time)) (t_grad) | $(frmt(val_lkl_time)) (t_val)"

			if isnan(val_loss) || isnan(batch_loss)
				error("Encountered invalid values in loss function.")
			end

			push!(history, :validation_likelihood, i, val_loss)
			push!(history, :val_lkl_time, i, val_lkl_time)

			if val_loss < best_val_loss
				best_val_loss = val_loss
				_patience = patience

				model = deepcopy(tr_model)
			else # else stop if the model has not improved for `patience` iterations
				_patience -= 1
				if _patience == 0
					@info "Stopped training after $(i) iterations, $((time() - start_time)/3600) hours."
					break
				end
			end
		end

		if time() - start_time > max_train_time # stop early if time is running out
			model = deepcopy(tr_model)
			@info "Stopped training after $(i) iterations, $((time() - start_time)/3600) hours due to time constraints."
			break
		end

		i += 1
	end

	(history=history, iterations=i, model=model, npars=sum(map(p -> length(p), ps)))
end


text:
export eigen2x2
function eigen2x2(A)
	Tp = eltype(A)
	a = A[1, 1]; b = A[1, 2]; c = A[2, 1]; d = A[2, 2];
	T = tr(A); D = det(A);

	L₁ = T / 2 + sqrt((T^2 / 4) - D)
	L₂ = T / 2 - sqrt((T^2 / 4) - D)
	values = [L₁, L₂]

	global vectors = []

	if !iszero(c)
		push!(vectors, [L₁ - d, c])
		push!(vectors, [L₂ - d, c])
	elseif !iszero(b)
		push!(vectors, [b, L₁ - a])
		push!(vectors, [b, L₂ - a])
	elseif iszero(b) && iszero(c)
		push!(vectors, [one(Tp), zero(Tp)])
		push!(vectors, [zero(Tp), one(Tp)])
	end

	return values, vectors
end


text:
# implementation of the GPUCompiler interfaces for generating SPIR-V code

# https://github.com/llvm/llvm-project/blob/master/clang/lib/Basic/Targets/SPIR.h
# https://github.com/KhronosGroup/LLVM-SPIRV-Backend/blob/master/llvm/docs/SPIR-V-Backend.rst
# https://github.com/KhronosGroup/SPIRV-LLVM-Translator/blob/master/docs/SPIRVRepresentationInLLVM.rst

const SPIRV_LLVM_Translator_jll = LazyModule("SPIRV_LLVM_Translator_jll", UUID("4a5d46fc-d8cf-5151-a261-86b458210efb"))
const SPIRV_Tools_jll = LazyModule("SPIRV_Tools_jll", UUID("6ac6d60f-d740-5983-97d7-a4482c0689f4"))

## target

export SPIRVCompilerTarget

Base.@kwdef struct SPIRVCompilerTarget <: AbstractCompilerTarget
end

llvm_triple(::SPIRVCompilerTarget) = Int===Int64 ? "spir64-unknown-unknown" : "spirv-unknown-unknown"

# SPIRV is not supported by our LLVM builds, so we can't get a target machine
llvm_machine(::SPIRVCompilerTarget) = nothing

llvm_datalayout(::SPIRVCompilerTarget) = Int===Int64 ?
    "e-i64:64-v16:16-v24:32-v32:32-v48:64-v96:128-v192:256-v256:256-v512:512-v1024:1024" :
    "e-p:32:32-i64:64-v16:16-v24:32-v32:32-v48:64-v96:128-v192:256-v256:256-v512:512-v1024:1024"

## job

# TODO: encode debug build or not in the compiler job
#       https://github.com/JuliaGPU/CUDAnative.jl/issues/368
runtime_slug(job::CompilerJob{SPIRVCompilerTarget}) = "spirv"

function process_kernel!(job::CompilerJob{SPIRVCompilerTarget}, mod::LLVM.Module, kernel::LLVM.Function)
    # calling convention
    for fun in functions(mod)
        callconv!(kernel, LLVM.API.LLVMSPIRFUNCCallConv)
    end
    callconv!(kernel, LLVM.API.LLVMSPIRKERNELCallConv)

    return kernel
end

function add_lowering_passes!(job::CompilerJob{SPIRVCompilerTarget}, pm::LLVM.PassManager)
    add!(pm, ModulePass("RemoveTrap", rm_trap!))
end

function mcgen(job::CompilerJob{SPIRVCompilerTarget}, mod::LLVM.Module, f::LLVM.Function,
               format=LLVM.API.LLVMAssemblyFile)
    # write the bitcode to a temporary file (the SPIRV Translator library doesn't have a C API)
    mktemp() do input, input_io
        write(input_io, mod)
        flush(input_io)

        # compile to SPIR-V
        mktemp() do output, output_io
            SPIRV_LLVM_Translator_jll.llvm_spirv() do translator
                cmd = `$translator`
                if format == LLVM.API.LLVMAssemblyFile
                    cmd = `$cmd -spirv-text`
                end
                cmd = `$cmd -o $output $input`
                run(cmd)
            end

            # read back the file
            if format == LLVM.API.LLVMAssemblyFile
                read(output_io, String)
            else
                read(output_io)
            end
        end
    end
end

# reimplementation that uses `spirv-dis`, giving much more pleasant output
function code_native(io::IO, job::CompilerJob{SPIRVCompilerTarget}; raw::Bool=false)
    obj, _ = codegen(:obj, job; strip=!raw)
    mktemp() do input_path, input_io
        write(input_io, obj)
        flush(input_io)

        SPIRV_Tools_jll.spirv_dis() do disassembler
            if io == stdout
                run(`$disassembler $input_path`)
            else
                mktemp() do output_path, output_io
                    run(`$disassembler $input_path -o $output_path`)
                    asm = read(output_io, String)
                    print(io, asm)
                end
            end
        end
    end
end

## LLVM passes

# SPIR-V does not support trap, and has no mechanism to abort compute kernels
# (OpKill is only available in fragment execution mode)
function rm_trap!(mod::LLVM.Module)
    job = current_job::CompilerJob
    changed = false
    @timeit_debug to "hide trap" begin

    if haskey(functions(mod), "llvm.trap")
        trap = functions(mod)["llvm.trap"]

        for use in uses(trap)
            val = user(use)
            if isa(val, LLVM.CallInst)
                unsafe_delete!(LLVM.parent(val), val)
                changed = true
            end
        end
    end

    end
    return changed
end


text:
using Quiqbox
using Quiqbox.Molden

mols = [
        ["H", "H"],
        ["N", "H", "H", "H"]
        ]
molNames = [
            "H2",
            "NH3"
            ]
br = 0.529177210903
# Data from CCCBDB: https://cccbdb.nist.gov
molCoords = [
             [[0.3705,0.0,0.0], [-0.3705,0.0,0.0]],
             [[0.0, 0.0, 0.1111], [0.0, 0.9316, -0.2592], [0.8068, -0.4658, -0.2592], [-0.8068, -0.4658, -0.2592]]
             ] ./ br

bfCoords = [molCoords..., GridBox(1, 1.2) |> gridCoords]
bfs = ["STO-3G"]
bsNames = push!(("-" .*molNames), "-Grid")
prefix = "Example"
for (nuc, nucCoords, molName, iMol) in zip(mols, molCoords, molNames, 1:length(mols)),
    (bfCoord, bsName) in zip(bfCoords[iMol:end], bsNames[iMol:end]),
    bf in bfs

    flag = (bfCoord == nucCoords)
    if flag
        nucConfig = [(bf, i) for i in nuc]
        bs = genBasisFunc.(bfCoord, nucConfig) |> flatten
    else
        bs = genBasisFunc.(bfCoord, bf) |> flatten
        bsName = "-Float"*bsName
    end

    # Number of spin-orbitals must not be smaller than numbers of electrons.
    fVars = try runHF(bs, nuc, nucCoords; printInfo=false) catch; continue end

    mol = Molecule(bs, nuc, nucCoords, fVars)
    fn = makeMoldenFile(mol; recordUMO=true, fileName=prefix*"_"*molName*"_"*bf*bsName)
end


text:
## Day 20: Trench Map ##
########################

day20_input = (10, 3)
day20_test = (4, 8)

# Given two player positions, simulates the dirac dice game
function simulate(a, b)
    dice = 10
    score_a = 0
    score_b = 0
    rolls = 0
    while true
        for _ = 1:3
            rolls += 1
            dice = wrap(dice + 1)
            a = wrap(a + dice)
        end
        score_a += a
        if score_a >= 1000
            return score_b, rolls
        end
        # @show score_a
        for _ = 1:3
            rolls += 1
            dice = wrap(dice + 1)
            b = wrap(b + dice)
        end
        score_b += b
        if score_b >= 1000
            return score_a, rolls
        end
        # @show score_b
    end
end

wrap(position) = position > 10 ? position - 10 : position

# Part 2, this looks more tricky. I thing the solution is to group universes
# that are identical and reduce them to a number.

# What is the state of the universe? (a, b, score_a, score_b). There are
# 10 * 10 * 31 * 31 = 96100 possible states (including 20 + 10 = 30 and 0 score)

function dirac(a, b)
    # To make indexing easier, players start with score 1 and win with 22 points.
    states = zeros(Int, (10, 10, 31, 31))
    states[a, b, 1, 1] = 1
    a_wins = 0
    b_wins = 0
    while sum(states) > 0
        @show "loop", sum(states), a_wins, b_wins
        # Evolve the metaverse
        new_states = zeros(Int, (10, 10, 31, 31))
        for a = 1:10, b = 1:10, score_a = 1:21, score_b = 1:21
            if states[a, b, score_a, score_b] > 0
                new_states .+= dirac_one_a(a, b, score_a, score_b) .* states[a, b, score_a, score_b]
            end
        end
        states = new_states
        # Filter out states that are already won and sum them up.
        for a = 1:10, b = 1:10, score_a = 22:31, score_b = 1:31
            a_wins += states[a, b, score_a, score_b]
            states[a, b, score_a, score_b] = 0
        end
        # Evolve the metaverse
        new_states = zeros(Int, (10, 10, 31, 31))
        for a = 1:10, b = 1:10, score_a = 1:21, score_b = 1:21
            if states[a, b, score_a, score_b] > 0
                new_states .+= dirac_one_b(a, b, score_a, score_b) .* states[a, b, score_a, score_b]
            end
        end
        states = new_states
        # Filter out states that are already won and sum them up.
        for a = 1:10, b = 1:10, score_a = 1:31, score_b = 22:31
            b_wins += states[a, b, score_a, score_b]
            states[a, b, score_a, score_b] = 0
        end
    end
    (a_wins, b_wins)
end

# You can roll 3*3*3 = 27 different universes which collapse into
roll_count = Dict(3 => 1, 4 => 3, 5 => 6, 6 => 7, 7 => 6, 8 => 3, 9 => 1)
@assert sum(values(roll_count)) == 27

function dirac_one_a(a, b, score_a, score_b)::Array
    states = zeros(Int, (10, 10, 31, 31))
    for roll in keys(roll_count)
        new_a = wrap(a + roll)
        states[new_a, b, score_a+new_a, score_b] = roll_count[roll]
    end
    states
end

function dirac_one_b(a, b, score_a, score_b)::Array
    states = zeros(Int, (10, 10, 31, 31))
    for roll in keys(roll_count)
        new_b = wrap(b + roll)
        states[a, new_b, score_a, score_b+new_b] = roll_count[roll]
    end
    states
end


text:
using OrdinaryDiffEq, DiffEqParamEstim, Distributions, Test,RecursiveArrayTools
      # ,BlackBoxOptim,

pf_func = function (du,u,p,t)
  du[1] = p[1] * u[1] - p[2] * u[1]*u[2]
  du[2] = -3.0 * u[2] + u[1]*u[2]
end
u0 = [1.0;1.0]
tspan = (0.0,10.0)
p = [1.5,1.0]
prob1 = ODEProblem(pf_func,u0,tspan,p)
sol = solve(prob1,Tsit5())
t = collect(range(0, stop=10, length=200))
function generate_data(sol,t)
  randomized = VectorOfArray([(sol(t[i]) + .01randn(2)) for i in 1:length(t)])
  data = convert(Array,randomized)
end
aggregate_data = convert(Array,VectorOfArray([generate_data(sol,t) for i in 1:100]))

distributions = [fit_mle(Normal,aggregate_data[i,j,:]) for i in 1:2, j in 1:200]
obj = build_loss_objective(prob1,Tsit5(),LogLikeLoss(t,distributions),
                                     maxiters=10000,verbose=false)
bound1 = Tuple{Float64, Float64}[(0.5, 5),(0.5, 5)]
result = bboptimize(obj;SearchRange = bound1, MaxSteps = 11e3)
@test result.archive_output.best_candidate ≈ [1.5,1.0] atol = 1e-1

data_distributions = [fit_mle(Normal,aggregate_data[i,j,:]) for i in 1:2, j in 1:200]
diff_distributions = [fit_mle(Normal,aggregate_data[i,j,:]-aggregate_data[i,j-1,:]) for j in 2:200, i in 1:2 ]
obj = build_loss_objective(prob1,Tsit5(),LogLikeLoss(t,data_distributions,diff_distributions),
                                     maxiters=10000,verbose=false)
bound1 = Tuple{Float64, Float64}[(0.5, 5),(0.5, 5)]
 result = bboptimize(obj;SearchRange = bound1, MaxSteps = 11e3)
@test result.archive_output.best_candidate ≈ [1.5,1.0] atol = 1e-1

data_distributions = [fit_mle(Normal,aggregate_data[i,j,:]) for i in 1:2, j in 1:200]
diff_distributions = [fit_mle(Normal,aggregate_data[i,j,:]-aggregate_data[i,j-1,:]) for j in 2:200, i in 1:2 ]
obj = build_loss_objective(prob1,Tsit5(),LogLikeLoss(t,data_distributions,diff_distributions,0.3),
                                     maxiters=10000,verbose=false)
bound1 = Tuple{Float64, Float64}[(0.5, 5),(0.5, 5)]
result = bboptimize(obj;SearchRange = bound1, MaxSteps = 11e3)
@test result.archive_output.best_candidate ≈ [1.5,1.0] atol = 1e-1

distributions = [fit_mle(MvNormal,aggregate_data[:,j,:]) for j in 1:200]
diff_distributions = [fit_mle(MvNormal,aggregate_data[:,j,:]-aggregate_data[:,j-1,:]) for j in 2:200]
priors = [Truncated(Normal(1.5,0.1),0,2),Truncated(Normal(1.0,0.1),0,1.5)]
obj = build_loss_objective(prob1,Tsit5(),LogLikeLoss(t,distributions,diff_distributions),
                                     maxiters=10000,verbose=false,priors=priors)
bound1 = Tuple{Float64, Float64}[(0.5, 5),(0.5, 5)]
result = bboptimize(obj;SearchRange = bound1, MaxSteps = 11e3)
@test result.archive_output.best_candidate ≈ [1.5,1.0] atol = 1e-1


text:
export nodetype, parentlinks

import AbstractTrees: nodetype, ParentLinks, StoredSiblings

function nodetype(::T) where {T <: AbstractJLBoostTree}
    T
end

function parentlinks(::T) where {T <: AbstractJLBoostTree}
    AbstractTrees.StoredParents()
end


text:
# Utility function for applying LayerNorm before a block
prenorm(planes, fn) = Chain(LayerNorm(planes), fn)

"""
    ChannelLayerNorm(sz::Int, λ = identity; ϵ = 1f-5)

A variant of LayerNorm where the input is normalised along the
channel dimension. The input is expected to have channel dimension with size
`sz`. It also applies a learnable shift and rescaling after the normalization.

Note that this is specifically for inputs with 4 dimensions in the format
(H, W, C, N) where H, W are the height and width of the input, C is the number
of channels, and N is the batch size.
"""
struct ChannelLayerNorm{F,D,T}
  λ::F
  diag::D
  ϵ::T
end

@functor ChannelLayerNorm

(m::ChannelLayerNorm)(x) = m.λ.(m.diag(MLUtils.normalise(x, dims = ndims(x) - 1, ϵ = m.ϵ)))

function ChannelLayerNorm(sz::Int, λ = identity; ϵ = 1f-5)
  diag = Flux.Diagonal(1, 1, sz)
  return ChannelLayerNorm(λ, diag, ϵ)
end


text:
include("simulate.jl")

test_parameters = Dict(
    "TrapFrequencies.axial_frequency" => 900e3,
    "TrapFrequencies.radial_frequency_1" => 1.3e6,
    "TrapFrequencies.radial_frequency_2" => 1.6e6,
)

num_ions = 1
b_field = 3.799953549
carrier_frequency = 2127405.1614560783

println("Testing a carrier Rabi flop:")
carrier_results = []
for duration in 0:1e-6:20e-6
    test_pulses = [
        Dict(
            "amp" => 1.0,
            "att" => 6.0,
            "dds_name" => "729G",
            "freq" => carrier_frequency,
            "phase" => 0.0,
            "time_off" => duration,
            "time_on" => 0.0
        ),
    ]

    result = simulate_with_ion_sim(test_parameters, test_pulses, num_ions, b_field)
    push!(carrier_results, result["D"])
end

println("Testing a blue sideband Rabi flop:")
sideband_results = []
for duration in 0:2e-6:50e-6
    test_pulses = [
        Dict(
            "amp" => 1.0,
            "att" => 6.0,
            "dds_name" => "729G",
            "freq" => carrier_frequency + test_parameters["TrapFrequencies.axial_frequency"],
            "phase" => 0.0,
            "time_off" => duration,
            "time_on" => 0.0
        ),
    ]

    result = simulate_with_ion_sim(test_parameters, test_pulses, num_ions, b_field)
    push!(sideband_results, result["D"])
end

println("Carrier results: $carrier_results")
println("Sideband results: $sideband_results")

max_stars_repo_name: rmshaffer/pulse-sequence-simulator

max_issues_repo_name: rmshaffer/pulse-sequence-simulator

max_forks_repo_name: rmshaffer/pulse-sequence-simulator


text:
module LeastSquaresSVM

using LinearAlgebra
using KernelFunctions
using Krylov
using Tullio
import MLJModelInterface

const MMI = MLJModelInterface
const BL = LinearAlgebra.BLAS

export SVM, LSSVC, LSSVR, KernelRBF, svmtrain, svmtrain_mc, svmpredict, LSSVClassifier,
LSSVRegressor

include("types.jl")
include("utils.jl")
include("tools.jl")
include("training.jl")
include("mlj_interface.jl")
# Some precompilation for the Tullio tools
include("precompile.jl")
_precompile_tullio()

end

max_stars_repo_name: edwinb-ai/LeastSquaresSVM

max_issues_repo_name: edwinb-ai/Elysivm

max_forks_repo_name: edwinb-ai/LeastSquaresSVM


text:
export SpikedWigner

struct SpikedWigner <: ContinuousMatrixDistribution
    beta::Integer
    n::Integer
    spikes::Array{Float64, 1}
    scaled::Bool
end

SpikedWigner(beta, n, spikes; scaled=false) = SpikedWigner(beta, n, spikes, scaled)
SpikedWigner(beta, n; scaled=false) = SpikedWigner(beta, n, [], scaled)

# SAMPLERS

function scaling(d::SpikedWigner)
    d.scaled ? 1/sqrt(d.n) : 1
end

function randreduced(d::SpikedWigner)
    if length(d.spikes) <= 1
        return randtridiagonal(d)
    else
        throw(">1 spike not yet implemented for Wigner matrices")
    end
end

function randtridiagonal(d::SpikedWigner)

    # diagonal and superdiagonal of B
    dv = rand(Normal(0, sqrt(2)), d.n)/sqrt(d.beta)
    ev = [rand(Chi(d.beta*(d.n-k))) for k in 1:(d.n-1)]/sqrt(d.beta)

    if length(d.spikes) == 1
        dv[1] += d.spikes[1] * sqrt(d.n)
    end

    SymTridiagonal(dv, ev) * scaling(d)
end


text:
struct Line
    start::CartesianIndex
    stop::CartesianIndex
end

function length(l::Line)
    return norm(l.stop - l.start)
end

function isvertical(l::Line)
    return l.start[1] == l.stop[1]
end

function ishorizontal(l::Line)
    return l.start[2] == l.stop[2]
end

function isdiagonal(l::Line)
    return abs(l.stop[1] - l.start[1]) == abs(l.stop[2] - l.start[2])
end

function getPointsOnLine(l::Line)
    points = CartesianIndex[]
    if isvertical(l)
        ystart = line.start[2] + 1
        ystop = line.stop[2] + 1
        step = 1
        if ystart > ystop
            step = -1
        end
        x = line.start[1] + 1
        y = ystart
        while y <= ystop
            points = push!(points, CartesianIndex(x,y))
            y += step
        end
    elseif ishorizontal(l)
        xstart = line.start[1] + 1
        xstop = line.stop[1] + 1
        step = 1
        if xstart > xstop
            step = -1
        end
        y = line.start[2] + 1
        x = xstart
        while x <= xstop
            points = push!(points, CartesianIndex(x,y))
            x += step
        end
    elseif isdiagonal(l)
        xstep = 1
        ystep = 1
        xstart = line.start[1] + 1
        ystart = line.start[2] + 1
        xstop = line.stop[1] + 1
        ystop = line.stop[2] + 1
        if xstop > xstart
            xstep = -1
        end
        if ystop > ystart
            ystep = -1
        end
        x = xstart
        y = ystart
        step = 1
        while step <= abs(xstop - xstart) + 1
            points = push!(points, CartesianIndex(x,y))
            x += xstep
            y += ystep
            step += 1
        end
    end
end

function solve(filename)
    #load txt file
    values = readlines(filename*".txt")
    lines = Line[]
    max_x = 0
    max_y = 0
    for line in values
        parts = split(line, '-')
        first = split(parts[1], ',')
        second = split(parts[2], ',')
        x1 = parse(Int, first[1])
        if x1 > max_x
            max_x = x1
        end
        x2string = replace(second[1], '>' => ' ')
        x2 = parse(Int, x2string)
        if x2 > max_x
            max_x = x2
        end
        y1 = parse(Int, first[2])
        if y1 > max_y
            max_y = y1
        end
        y2 = parse(Int, second[2])
        if y2 > max_y
            max_y = y2
        end
        start = CartesianIndex(x1, y1)
        stop = CartesianIndex(x2, y2)
        lines = push!(lines, Line(start, stop))
    end

    map = zeros(max_x + 1, max_y + 1)
    for line in lines
        #println(line.start)
        #println(line.stop)
        if isvertical(line)
            #println("vertical")
            ystart = line.start[2] + 1
            ystop = line.stop[2] + 1
            if ystart > ystop
                tmp = ystart
                ystart = ystop
                ystop = tmp
            end
            x = line.start[1] + 1
            y = ystart
            while y <= ystop
                map[x,y] += 1
                y += 1
            end
        elseif ishorizontal(line)
            #println("horizontal")
            xstart = line.start[1] + 1
            xstop = line.stop[1] + 1
            if xstart > xstop
                tmp = xstart
                xstart = xstop
                xstop = tmp
            end
            y = line.start[2] + 1
            x = xstart
            while x <= xstop
                #println("mark: ("*string(x)*", "*string(y)*")")
                map[x,y] += 1
                #println(map[x,y])
                x += 1
            end
        end
    end
    count_morethantwo = 0
    for el in map
        if el >= 2
            count_morethantwo += 1
        end
    end
    println(count_morethantwo)
end

solve("smalltest")
solve("input")

function solve2(filename)
    #load txt file
    values = readlines(filename*".txt")
    lines = Line[]
    max_x = 0
    max_y = 0
    for line in values
        parts = split(line, '-')
        first = split(parts[1], ',')
        second = split(parts[2], ',')
        x1 = parse(Int, first[1])
        if x1 > max_x
            max_x = x1
        end
        x2string = replace(second[1], '>' => ' ')
        x2 = parse(Int, x2string)
        if x2 > max_x
            max_x = x2
        end
        y1 = parse(Int, first[2])
        if y1 > max_y
            max_y = y1
        end
        y2 = parse(Int, second[2])
        if y2 > max_y
            max_y = y2
        end
        start = CartesianIndex(x1, y1)
        stop = CartesianIndex(x2, y2)
        lines = push!(lines, Line(start, stop))
    end

    map = zeros(max_x + 1, max_y + 1)
    for line in lines
        if isvertical(line)
            ystart = line.start[2] + 1
            ystop = line.stop[2] + 1
            if ystart > ystop
                tmp = ystart
                ystart = ystop
                ystop = tmp
            end
            x = line.start[1] + 1
            y = ystart
            while y <= ystop
                map[x,y] += 1
                y += 1
            end
        elseif ishorizontal(line)
            xstart = line.start[1] + 1
            xstop = line.stop[1] + 1
            if xstart > xstop
                tmp = xstart
                xstart = xstop
                xstop = tmp
            end
            y = line.start[2] + 1
            x = xstart
            while x <= xstop
                map[x,y] += 1
                x += 1
            end
        elseif isdiagonal(line)
            x_right = false
            y_down = false
            xstart = line.start[1] + 1
            ystart = line.start[2] + 1
            xstop = line.stop[1] + 1
            ystop = line.stop[2] + 1
            if xstop > xstart
                x_right = true
            end
            if ystop > ystart
                y_down = true
            end
            xstep = 1
            ystep = 1
            if x_right == false
                xstep = -1
            end
            if y_down ==false
                ystep = -1
            end
            x = xstart
            y = ystart
            step = 1
            while step <= abs(xstop - xstart) + 1
                map[x,y] += 1
                x += xstep
                y += ystep
                step += 1
            end
        end
    end
    count_morethantwo = 0
    for el in map
        if el >= 2
            count_morethantwo += 1
        end
    end
    println(count_morethantwo)
end

solve2("smalltest")
solve2("input")


text:
include("concordance.jl")

using DataFrames, JuMP, Ipopt, DelimitedFiles;

#Change GFCF_By_Industry_Asset to 19 sector
flows97 = ExcelReaders.readxlsheet("data"*pathmark*"flow1997.xls", "180x22Combined");
flows = DataFrame(flows97[4:182, 4:25], :auto);
dataTypeEx = flows[41,4]
for i in [1:1:179;]
    for j in [1:1:22;]
        if typeof(flows[i,j]) == typeof(dataTypeEx)
            flows[i,j] = Float64(0.0)
        end
    end
end

# Aggregating to 19 sector
flows.x5 = flows.x5+flows.x6+flows.x7;
flows = select!(flows, Not(:x6));
flows = select!(flows, Not(:x7));
flows.x10 = flows.x10+flows.x11;
flows = select!(flows, Not(:x11));
flows.x15 = flows.x15+flows.x16;
flows = select!(flows, Not(:x16));

# Renaming to 19 Sectors codes
rename!(flows, :x1 => :A)
rename!(flows, :x2 => :B)
rename!(flows, :x3 => :D)
rename!(flows, :x4 => :E)
rename!(flows, :x5 => :C)
rename!(flows, :x8 => :F)
rename!(flows, :x9 => :G)
rename!(flows, :x10 => :I)
rename!(flows, :x12 => :J)
rename!(flows, :x13 => :K)
rename!(flows, :x14 => :L)
rename!(flows, :x15 => :M)
rename!(flows, :x17 => :N)
rename!(flows, :x18 => :P)
rename!(flows, :x19 => :Q)
rename!(flows, :x20 => :R)
rename!(flows, :x21 => :H)
rename!(flows, :x22 => :S)

rowCode = flows97[4:182,2];
containsSpace = findall( x -> occursin(" ", x), rowCode);
rowCode[containsSpace] = replace.(rowCode[containsSpace], " " => "");
rowCode19 = Array{Union{Nothing, String}}(nothing, length(rowCode));
for i in eachindex(rowCode);
    rowCode19[i] = Comm180To19[rowCode[i]];
end
insertcols!(flows ,1, :Industry => rowCode19);
insertcols!(flows ,15, :O => zeros(179));
splitIndustry = groupby(flows, :Industry);
flows = combine(splitIndustry, valuecols(splitIndustry) .=> sum);

push!(flows, ["A" zeros(1,19)])
push!(flows, ["D" zeros(1,19)])
push!(flows, ["H" zeros(1,19)])
push!(flows, ["K" zeros(1,19)])
push!(flows, ["N" zeros(1,19)])
push!(flows, ["P" zeros(1,19)])
push!(flows, ["Q" zeros(1,19)])
push!(flows, ["R" zeros(1,19)])
push!(flows, ["S" zeros(1,19)])

sort!(flows)

flows = select!(flows, Not(:Industry));
insertcols!(flows ,20, :Dwelling => zeros(19));
push!(flows, zeros(1,20));

ausGFCF = ExcelReaders.readxlsheet("data"*pathmark*"5204064_GFCF_By_Industry_Asset.xls", "Data1");
ausGFCF = ausGFCF[:,Not(findall(x -> occursin("ALL INDUSTRIES ;", x), string.(ausGFCF[1,:])))]
ausGFCF2019 = ausGFCF[findall(x -> occursin("2019", x), string.(ausGFCF[:,1])),
    findall(x -> occursin("Gross fixed capital formation: Current prices ;", x), string.(ausGFCF[1,:]))];
ausGFCF2019=[ausGFCF2019 (ausGFCF[findall(x -> occursin("2019", x), string.(ausGFCF[:,1])),
findall(x -> occursin("Dwellings: Current prices ;", x), string.(ausGFCF[1,:]))]
+ausGFCF[findall(x -> occursin("2019", x), string.(ausGFCF[:,1])),
findall(x -> occursin("Ownership transfer costs: Current prices ;", x), string.(ausGFCF[1,:]))])];

IOIGAs19=Array{Union{Nothing, String}}(nothing, length(IOIG));
for i in eachindex(IOIG);
    IOIGAs19[i] = IOIGTo19[IOIG[i]]
end

IOSource8 = ExcelReaders.readxlsheet("data"*pathmark*"5209055001DO001_201819.xls", "Table 8");

privateGFCF = IOSource8[4:117,findall(x -> occursin("Private ; Gross Fixed Capital Formation", x), string.(IOSource8[2,:]))]
publicGFCF = IOSource8[4:117,findall(x -> occursin("Public Corporations ; Gross Fixed Capital Formation", x), string.(IOSource8[2,:]))]
govGFCF = IOSource8[4:117,findall(x -> occursin("General Government ; Gross Fixed Capital Formation", x), string.(IOSource8[2,:]))]

privateGFCFTot = IOSource8[127,findall(x -> occursin("Private ; Gross Fixed Capital Formation", x), string.(IOSource8[2,:]))]
publicGFCFTot = IOSource8[127,findall(x -> occursin("Public Corporations ; Gross Fixed Capital Formation", x), string.(IOSource8[2,:]))]
govGFCFTot = IOSource8[127,findall(x -> occursin("General Government ; Gross Fixed Capital Formation", x), string.(IOSource8[2,:]))]

privateGFCFSum = sum(privateGFCF)
publicGFCFSum = sum(publicGFCF)
govGFCFSum = sum(govGFCF)

for i in eachindex(privateGFCF);
    privateGFCF[i] = privateGFCF[i]*privateGFCFTot/privateGFCFSum
    publicGFCF[i] = publicGFCF[i]*publicGFCFTot/publicGFCFSum
    govGFCF[i] = govGFCF[i]*govGFCFTot/govGFCFSum
end

ausCapReceiv = DataFrame(privateGFCF+publicGFCF+govGFCF, [:GFCF]);
insertcols!(ausCapReceiv ,1, :Industry => IOIGAs19);
splitIndustry = groupby(ausCapReceiv, :Industry);
ausCapReceiv = combine(splitIndustry, valuecols(splitIndustry) .=> sum);
sort!(ausCapReceiv);

govCapReceiv = DataFrame(govGFCF, [:GFCF]);
insertcols!(govCapReceiv ,1, :Industry => IOIGAs19);
splitIndustry = groupby(govCapReceiv, :Industry);
govCapReceiv = combine(splitIndustry, valuecols(splitIndustry) .=> sum);
sort!(govCapReceiv);

ausCapReceivable = [ausCapReceiv.GFCF_sum; 0];
govCapReceivable = [govCapReceiv.GFCF_sum; 0];

#balancing rows and collumn sums(very small diff)
ausCapReceivableSum = sum(ausCapReceivable)
ausGFCF2019Sum = sum(ausGFCF2019)
avSum = (ausCapReceivableSum + ausGFCF2019Sum)/2;
for i in eachindex(ausCapReceivable)
    ausCapReceivable[i]=ausCapReceivable[i]*avSum/ausCapReceivableSum;
end
for i in eachindex(ausGFCF2019)
    ausGFCF2019[i]=ausGFCF2019[i]*avSum/ausGFCF2019Sum;
end

modCap = Model(Ipopt.Optimizer);
@variable(modCap, x[1:length(ausCapReceivable), 1:length(ausGFCF2019)]);
@NLobjective(modCap, Min, sum(((x[i,j] - (flows[i,j]+1000))/ (flows[i,j]+1000))^ 2 for i in eachindex(ausCapReceivable), j in eachindex(ausGFCF2019)));
for j in eachindex(ausGFCF2019)
    @constraint(modCap, sum(x[:,j]) == ausGFCF2019[j]+20000);
end;
for k in eachindex(ausCapReceivable)
    @constraint(modCap, sum(x[k,:]) == ausCapReceivable[k]+20000);
end;
optimize!(modCap);

capitalFlowsTitleRow = permutedims([ANZSICDivShort; "Dwelling"]);
capitalFlowsTitleCol = ["Titles"; ANZSICDivShort; "Dwelling"];
capitalFlowsDraft =[capitalFlowsTitleCol [capitalFlowsTitleRow; value.(x)]]

writedlm("data"*pathmark*"capitalFlowsDraft.csv", ',');

max_stars_repo_name: uq-aibe/Capital

max_issues_repo_name: uq-aibe/Capital

max_forks_repo_name: uq-aibe/Capital


text:
module MultivariateSeries

using LinearAlgebra
using MultivariatePolynomials
using DynamicPolynomials

include("polynomials.jl")
include("series.jl")
include("moments.jl")
include("hankel.jl")
include("newton.jl")
include("diagonalisation.jl")
include("decompose.jl")

export Seq, seq

mutable struct Seq{T} val::T end

function seq(args...)
    if length(args)>1
        Seq([args...])
    else
        Seq(args[1])
    end
end

function DynamicPolynomials.MonomialVector(V::Vector{PolyVar{true}}, rg::Seq)
     L = DynamicPolynomials.Monomial{true}[]
     for i in rg.val
         append!(L, DynamicPolynomials.monomials(V,i))
     end
     L
end

function DynamicPolynomials.MonomialVector(V::Vector{PolyVar{true}}, rg::UnitRange{Int64})
     L = DynamicPolynomials.Monomial{true}[]
     for i in rg
         append!(L, DynamicPolynomials.monomials(V,i))
     end
     L
end

end


text:
using HDF5

@testset "Output tests:" begin

    umin   = 0.0
    umax   = 2.0
    unodes = 16

    ucoord = GaussLobatto("u", umin, umax, unodes)

    xmin   = -10.0
    xmax   =  10.0
    xnodes =  20

    xcoord = Cartesian{2}("x", xmin, xmax, xnodes, endpoint=false)

    ymin   = -20.0
    ymax   =  20.0
    ynodes =  20

    ycoord = Cartesian{3}("y", ymin, ymax, ynodes, endpoint=false)

    f     = [0.5 * x1.^2 .* cos.(x2) .* sin.(x3) for x1 in ucoord[:], x2 in xcoord[:], x3 in ycoord[:]]
    g     = [666 for x1 in ucoord[:], x2 in xcoord[:], x3 in ycoord[:]]

    chart  = Jecco.Chart(ucoord, xcoord, ycoord)

    mins      = Jecco.min(chart)
    deltas    = Jecco.delta(chart)
    maxs      = Jecco.max(chart)
    gridtypes = Jecco.coord_type(chart)
    names     = Jecco.name(chart)

    field1 = Jecco.Field("f", f, chart)
    field2 = Jecco.Field("g", g, chart)

    # write the contents
    tinfo  = Jecco.TimeInfo(1, 10.0, 0.1, 0.0)
    dir    = tempname()
    prefix = "data_"
    out    = Jecco.Output(dir, prefix, tinfo)
    out(field1, field2)

    # and now read the file back in
    fn  = "$(prefix)00000001.h5"

    fid = h5open(dir * "/" * fn, "r")

    grp = fid["data/1"]

    @test read_attribute(grp, "time") == 10.0
    @test read_attribute(grp, "dt")   == 0.1

    dset1 = grp["fields/f"]
    data1 = read(dset1)

    @test all(data1 .== f)
    @test all(read_attribute(dset1, "gridGlobalOffset") .== mins[end:-1:1])
    @test all(read_attribute(dset1, "gridMax")          .== maxs[end:-1:1])
    @test all(read_attribute(dset1, "gridType")         .== gridtypes[end:-1:1])
    @test all(read_attribute(dset1, "axisLabels")       .== names[end:-1:1])

    dset2 = grp["fields/g"]
    data2 = read(dset2)

    @test all(data2 .== g)

    close(fid)
end


text:
# ### Hisogram plot 2d

#md # ![FILE_NAME.png](images/FILE_NAME.png)

x = randn(50000);
y = x * 1.2 + randn(50000);

plot(
    x,
    y,
    seriestype = :histogram2d,
    c = :vik,
    nbins = 200,
    show_empty_bins = :true
)

## equivalent
## histogram2d(
##     x,
##     y,
##     c = :vik,
##     nbins = 200,
##     show_empty_bins = :true
## )


text:
mutable struct PsdConeTriangleLOBPCG{T} <: AbstractConvexCone{T}
    dim::Int
    sqrt_dim::Int
    iteration::Int
    is_subspace_positive::Bool
    X::Matrix{T}
    U::Matrix{T}
    lobpcg::LOBPCGIterable{T}
    exact_projections::Int
    lobpcg_iterations::Int
    max_iter::Int
    subspace_dim_history::Vector{Int}

    function PsdConeTriangleLOBPCG{T}(dim::Int; buffer_size::Int=-1, max_iter::Int=10, max_subspace_dim::Int=-1, verbosity::Int=0) where {T}
        dim >= 0       || throw(DomainError(dim, "dimension must be nonnegative"))
        sqrt_dim = Int(1 / 2 * (sqrt(8 * dim + 1) - 1)) # Solution of (sqrt_dim^2 + sqrt_dim)/2 = length(x) obtained by WolframAlpha
        sqrt_dim * (sqrt_dim + 1) / 2 == dim || throw(DomainError(dim, "dimension must be a square"))

        if buffer_size < 0
            buffer_size = Int(floor(min(max(sqrt_dim / 50, 3), 20)))
        end
        if max_subspace_dim < 0
            max_subspace_dim = Int(floor(sqrt_dim/4))
        end

        lobpcg = LOBPCGIterable(zeros(T, sqrt_dim, sqrt_dim), verbosity=verbosity,
            buffer_size=buffer_size, λ_barrier=T(0), max_dim=max_subspace_dim, min_full_iters=1)

        new(dim, sqrt_dim, 0, true,
            zeros(T, sqrt_dim, sqrt_dim), # X
            zeros(T, sqrt_dim, 0), # U
            lobpcg,
            0, 0, max_iter, # exact_projections, lobpcg_iterations
            zeros(Int, 0) # subspace_dim_history
        )
    end
end
PsdConeTriangleLOBPCG(args...; kwargs...) = PsdConeTriangleLOBPCG{DefaultFloat}(args...; kwargs...)

function get_tolerance(cone::PsdConeTriangleLOBPCG{T}) where T
    return T(max(sqrt(cone.sqrt_dim) / cone.iteration^(1.01) * 10, 1e-7))
end

function project_exact!(x::AbstractArray{T}, cone::PsdConeTriangleLOBPCG{T}) where {T}
    cone.exact_projections += 1
    populate_upper_triangle!(cone.X, x)
    λ, U  = eigen!(Symmetric(cone.X))
    zero_tol = 1e-9
    cone.is_subspace_positive = sum(λ .> zero_tol) < cone.sqrt_dim/2
    if cone.is_subspace_positive
        range = max(sum(λ .<= zero_tol) - cone.lobpcg.buffer_size + 1, 1):cone.sqrt_dim
    else
        range = 1:min(sum(λ .< -zero_tol) + cone.lobpcg.buffer_size, cone.sqrt_dim)
        λ .= -λ
        populate_upper_triangle!(cone.X, x) # Restore cone.X altered by eigen!
    end
    cone.U = view(U, :, range)
    construct_projection!(x, cone, view(λ, range))
end

function construct_projection!(x, cone::PsdConeTriangleLOBPCG{T}, λ) where {T}
    append!(cone.subspace_dim_history, size(cone.U, 2))
    scaled_U = rmul!(copy(cone.U), Diagonal(sqrt.(max.(λ, 0))))
    if cone.is_subspace_positive
        BLAS.syrk!('U', 'N', one(T), scaled_U, zero(T), cone.X)
    else
        BLAS.syrk!('U', 'N', one(T), scaled_U, one(T), cone.X)
    end
    extract_upper_triangle!(cone.X, x)
end

function project!(x::AbstractArray, cone::PsdConeTriangleLOBPCG{T}) where {T}
    sqrt_dim = cone.sqrt_dim
    cone.iteration += 1

    if size(cone.U, 2) <= cone.lobpcg.max_dim && cone.iteration > 1
        populate_upper_triangle!(cone.X, x)
        initialize!(cone.lobpcg, Symmetric(cone.X), cone.U, is_orthonormal=true,
            which = cone.is_subspace_positive ? :largest : :smallest)
        cone.lobpcg.tol = get_tolerance(cone)
        try
            cone.lobpcg, status = lobpcg!(cone.lobpcg, cone.max_iter)
            cone.lobpcg_iterations += cone.lobpcg.iteration - 1
        catch e
            if isa(e, PosDefException)
                status = :error
                @warn "LOBPCG failed; switching to exact eigendecompositon."
            else
                throw(e)
            end
        end
    else
        status = :not_called
    end

    if status != :converged
        return project_exact!(x, cone)
    else
        cone.U = cone.lobpcg.X
        if !cone.is_subspace_positive
            cone.lobpcg.λ .= -cone.lobpcg.λ
        end
        construct_projection!(x, cone, cone.lobpcg.λ)
    end
end

function reset_iteration_counters!(cone::PsdConeTriangleLOBPCG)
    cone.exact_projections = 0
    cone.lobpcg_iterations = 0
    cone.iteration = 1
end

allocate_memory!(cone::PsdConeTriangleLOBPCG) = nothing


text:
function Base.show(io::IO, ::MIME"text/plain", g::Gene)
    for fname in fieldnames(Gene)
        _pretty_print_keyvals(io, "Gene.$(string(fname)): ", getfield(g, fname))
    end
end


text:
module StaggeredGridCalculus

# @reexport makes all exported symbols of the exported packages available in module using
# StaggeredGridCalculus.
using Reexport
@reexport using LinearAlgebra, SparseArrays
using AbbreviatedTypes
using Base.Threads: @threads  # used in matrixfree/differential/*.jl

# The order of inclusion matters: if types or functions in file A are used in file B, file A
# must be included first.
include("enumtype.jl")
include("util.jl")
include("grid.jl")
# include("gridgen.jl")
include("matrix/matrix.jl")
include("matrixfree/matrixfree.jl")
include("pml.jl")

end # module StaggeredGridCalculus


text:
using Agents, Random

mutable struct SheepWolf <: AbstractAgent
    id::Int
    pos::Dims{2}
    type::Symbol # :sheep or :wolf
    energy::Float64
    reproduction_prob::Float64
    Δenergy::Float64
end

# Simple helper functions
Sheep(id, pos, energy, repr, Δe) = SheepWolf(id, pos, :sheep, energy, repr, Δe)
Wolf(id, pos, energy, repr, Δe) = SheepWolf(id, pos, :wolf, energy, repr, Δe)

"""
``` julia
predator_prey(;
    n_sheep = 100,
    n_wolves = 50,
    dims = (20, 20),
    regrowth_time = 30,
    Δenergy_sheep = 4,
    Δenergy_wolf = 20,
    sheep_reproduce = 0.04,
    wolf_reproduce = 0.05,
    seed = 23182,
)
```
Same as in [Predator-prey dynamics](@ref).
"""
function predator_prey(;
    n_sheep = 100,
    n_wolves = 50,
    dims = (20, 20),
    regrowth_time = 30,
    Δenergy_sheep = 4,
    Δenergy_wolf = 20,
    sheep_reproduce = 0.04,
    wolf_reproduce = 0.05,
    seed = 23182,
)

    rng = MersenneTwister(seed)
    space = GridSpace(dims, periodic = false)
    ## Model properties contain the grass as two arrays: whether it is fully grown
    ## and the time to regrow. Also have static parameter `regrowth_time`.
    ## Notice how the properties are a `NamedTuple` to ensure type stability.
    properties = (
        fully_grown = falses(dims),
        countdown = zeros(Int, dims),
        regrowth_time = regrowth_time,
    )
    model = ABM(SheepWolf, space; properties, rng, scheduler = random_activation)
    id = 0
    for _ in 1:n_sheep
        id += 1
        energy = rand(1:(Δenergy_sheep*2)) - 1
        sheep = Sheep(id, (0, 0), energy, sheep_reproduce, Δenergy_sheep)
        add_agent!(sheep, model)
    end
    for _ in 1:n_wolves
        id += 1
        energy = rand(1:(Δenergy_wolf*2)) - 1
        wolf = Wolf(id, (0, 0), energy, wolf_reproduce, Δenergy_wolf)
        add_agent!(wolf, model)
    end
    for p in positions(model) # random grass initial growth
        fully_grown = rand(model.rng, Bool)
        countdown = fully_grown ? regrowth_time : rand(model.rng, 1:regrowth_time) - 1
        model.countdown[p...] = countdown
        model.fully_grown[p...] = fully_grown
    end
    return model
end

function sheepwolf_step!(agent::SheepWolf, model)
    if agent.type == :sheep
        sheep_step!(agent, model)
    else # then `agent.type == :wolf`
        wolf_step!(agent, model)
    end
end

function sheep_step!(sheep, model)
    walk!(sheep, rand, model)
    sheep.energy -= 1
    sheep_eat!(sheep, model)
    if sheep.energy < 0
        kill_agent!(sheep, model)
        return
    end
    if rand(model.rng) <= sheep.reproduction_prob
        wolfsheep_reproduce!(sheep, model)
    end
end

function wolf_step!(wolf, model)
    walk!(wolf, rand, model)
    wolf.energy -= 1
    agents = collect(agents_in_position(wolf.pos, model))
    dinner = filter!(x -> x.type == :sheep, agents)
    wolf_eat!(wolf, dinner, model)
    if wolf.energy < 0
        kill_agent!(wolf, model)
        return
    end
    if rand(model.rng) <= wolf.reproduction_prob
        wolfsheep_reproduce!(wolf, model)
    end
end

function sheep_eat!(sheep, model)
    if model.fully_grown[sheep.pos...]
        sheep.energy += sheep.Δenergy
        model.fully_grown[sheep.pos...] = false
    end
end

function wolf_eat!(wolf, sheep, model)
    if !isempty(sheep)
        dinner = rand(model.rng, sheep)
        kill_agent!(dinner, model)
        wolf.energy += wolf.Δenergy
    end
end

function wolfsheep_reproduce!(agent, model)
    agent.energy /= 2
    id = nextid(model)
    offspring = SheepWolf(
        id,
        agent.pos,
        agent.type,
        agent.energy,
        agent.reproduction_prob,
        agent.Δenergy,
    )
    add_agent_pos!(offspring, model)
    return
end

function grass_step!(model)
    @inbounds for p in positions(model) # we don't have to enable bound checking
        if !(model.fully_grown[p...])
            if model.countdown[p...] ≤ 0
                model.fully_grown[p...] = true
                model.countdown[p...] = model.regrowth_time
            else
                model.countdown[p...] -= 1
            end
        end
    end
end


text:
abstract type Node end

for pubsub in ((:publisher, :pubs), (:subscriber, :subs))
    @eval $(Symbol("get", pubsub[1]))(node::Node, index::Integer) =
        getIO(node).$(pubsub[2])[index]
    @eval function $(Symbol("get", pubsub[1]))(node::Node, name::String)
        index = findfirst(getIO(node).$(pubsub[2])) do msg
            getname(msg) == name
        end
        if !isnothing(index)
            return $(Symbol("get", pubsub[1]))(node, index)
        end
        return nothing
    end
    @eval $(Symbol("num", pubsub[1], "s"))(node::Node) = length(getIO(node).$(pubsub[2]))
end

# %%
for pubsub in ((:publisher, :pubs), (:subscriber, :subs))
    println(pubsub)

    # @eval $(Symbol("get", pubsub[1]))(node::Node, index::Integer) =
    #     getIO(node).$(pubsub[2])[index]
    println((Symbol("get", pubsub[1])))
end


text:
# Julia wrapper for header: /opt/local/include/gdalexif.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function TIFFSwabShort(arg1::Ptr{Cint})
    ccall((:TIFFSwabShort,libgdal),Void,(Ptr{Cint},),arg1)
end

function TIFFSwabLong(arg1::Ptr{Cint})
    ccall((:TIFFSwabLong,libgdal),Void,(Ptr{Cint},),arg1)
end

function TIFFSwabDouble(arg1::Ptr{Cdouble})
    ccall((:TIFFSwabDouble,libgdal),Void,(Ptr{Cdouble},),arg1)
end

function TIFFSwabArrayOfShort(arg1::Ptr{Cint},arg2::Culong)
    ccall((:TIFFSwabArrayOfShort,libgdal),Void,(Ptr{Cint},Culong),arg1,arg2)
end

function TIFFSwabArrayOfLong(arg1::Ptr{Cint},arg2::Culong)
    ccall((:TIFFSwabArrayOfLong,libgdal),Void,(Ptr{Cint},Culong),arg1,arg2)
end

function TIFFSwabArrayOfDouble(arg1::Ptr{Cdouble},arg2::Culong)
    ccall((:TIFFSwabArrayOfDouble,libgdal),Void,(Ptr{Cdouble},Culong),arg1,arg2)
end


text:
abstract type AbstractParameter end
abstract type AbstractParameterSet end
abstract type ParameterSpace<:AbstractParameterSet end
abstract type ParameterGroup<:AbstractParameterSet end

struct Parameter{T,S<:String} <: AbstractParameter
    value::T
    description::S
end
struct BoundedParameter{T,S<:String} <: AbstractParameter
    value::T
    description::S
    lb::T
    ub::T
end

function Base.getproperty(name::ParameterGroup,val::Symbol)
    if typeof(getfield(name,val))<:AbstractParameter
        return getfield(getfield(name,val),:value)
    else
        return getfield(name,val)
    end
end

function description(name::ParameterGroup,val::Symbol)
    if typeof(getfield(name,val))<:AbstractParameter
        return getfield(getfield(name,val),:description)
    else
        @warn "No description provided"
        return nothing
    end
end

function bounds(name::ParameterGroup,val::Symbol)
    if typeof(getfield(name,val))<:BoundedParameter
        return getfield(getfield(name,val),:lb), getfield(getfield(name,val),:ub)
    else
        @warn "No bounds defined"
        return nothing
    end
end