text:
using Luna
import PyPlot: plt

#= In this example we simulate degenerate four-wave mixing between a circularly polarised
pump pulse at 400 nm and a seed pulse at 800 nm to generate a circularly polarised idler
at 266 nm. By switching the helicity of the seed relative to the pump, we can suppress
the idler generation.
The parameters are chosen to correspond to those in:
A. Lekosiotis, et al.,
Generation of broadband circularly polarized deep-ultraviolet pulses in hollow capillary fibers,
Opt. Lett., 45, 5648, (2020),
DOI: 10.1364/OL.400362.
=#

λp = 400e-9 # central wavelength of pump pulse
ep = 100e-6 # energy in pump pulse

λs = 800e-9 # central wavelength of seed pulse
es = 27e-6 # energy in seed pulse

τfwhm = 30e-15 # pulse duration of both pulses

radius = 75e-6 # HCF core radius
flength = 1.3 # HCF length
gas = :He
pressure = 1.8 # helium pressure in bar

pump = Pulses.GaussPulse(;λ0=λp, energy=ep, τfwhm, polarisation=:circular)
seed_equal = Pulses.GaussPulse(;λ0=λs, energy=es, τfwhm, polarisation=1.0)
seed_opposite = Pulses.GaussPulse(;λ0=λs, energy=es, τfwhm, polarisation=-1.0)

# Here we are not specifiying a mode, so Luna will automatically choose to propagate in
# the fundamental mode only--but because the pulses are circularly polarised, two modes
# with perpendicular linear polarisation are used.
equal = prop_capillary(radius, flength, gas, pressure; λ0=λp, pulses=[pump, seed_equal],
                       trange=500e-15, λlims=(200e-9, 1500e-9))

opposite = prop_capillary(radius, flength, gas, pressure; λ0=λp, pulses=[pump, seed_opposite],
                          trange=500e-15, λlims=(200e-9, 1500e-9))

# retrieve the spectral energy density at the output
λ, Iλequal = Processing.getIω(equal, :λ, flength)
_, Iλopposite = Processing.getIω(opposite, :λ, flength)

# sum over the two modes
Iλequal = dropdims(sum(Iλequal; dims=2); dims=(2, 3))
Iλopposite = dropdims(sum(Iλopposite; dims=2); dims=(2, 3))

# plot the result
plt.figure()
plt.plot(λ*1e9, Iλequal*1e-3, label="Equal helicity")
plt.plot(λ*1e9, Iλopposite*1e-3, label="Opposite helicity")
plt.xlim(200, 950)
plt.xlabel("Wavelength (nm)")
plt.ylabel("Spectral energy density (μJ/nm)")
plt.ylim(ymin=0)
plt.legend()


text:
function rectangles!(r)

    function change(side, num)
        return HorizonSide((Int(side) + num) % 4)
    end

    function go_around(side)
        num = 0
        while isborder(r, side)
            if !isborder(r, change(side, 1))
                move!(r, change(side, 1))
                num += 1
            else
                break
            end
        end
        if isborder(r, side)
            for i in 1:num
                if !isborder(r, change(side, 3))
                    move!(r, change(side, 3))
                end
            end
            return false
        else
            move!(r, side)
            while isborder(r, change(side, 3))
                move!(r, side)
            end
            for i in 1:num
                if !isborder(r, change(side, 3))
                    move!(r, change(side, 3))
                end
            end
            return true
        end
    end

    function move_and_kol(side, num)
        flag = false
        ok = true
        count = 0
        while true
            if (isborder(r, change(side, num)) && !flag)
                flag = true
                count += 1
            elseif !isborder(r, change(side, num))
                flag = false
            end
            if (isborder(r, side))
                ok = go_around(side)
            end
            if ok
                if !isborder(r, side)
                    move!(r, side)
                end
            else
                break
            end
        end
        return count
    end

    function solve()
        side = Ost
        n = 1
        count = 0
        while !isborder(r, Nord)
            count += move_and_kol(side, n)
            side = change(side, 2)
            n = (n + 2) % 4
            move!(r, Nord)
        end
        println(count)
    end

    solve()
end

ext: jl

max_stars_repo_name: dmitrii-a-piskunov/mirea-progs

max_issues_repo_name: dmitrii-a-piskunov/mirea-progs

max_forks_repo_name: dmitrii-a-piskunov/mirea-progs


text:
# =================================================================================================
# =                                         Hyland83a.jl                                          =
# =================================================================================================

#=
Testing the virial coefficients. Table 1 from reference [2].
=#
# Baa
@test -55.94 ≈ Psychro.Baa(173.15)*1e6 atol=0.01
@test -13.15 ≈ Psychro.Baa(273.15)*1e6 atol=0.01
@test 3.72 ≈ Psychro.Baa(373.15)*1e6 atol=0.01
@test 12.31 ≈ Psychro.Baa(473.15)*1e6 atol=0.01

# dBaa/dT
@test 0.7240 ≈ Psychro.dBaa(173.15)*1e6 atol=0.0001
@test 0.2460 ≈ Psychro.dBaa(273.15)*1e6 atol=0.0001
@test 0.1146 ≈ Psychro.dBaa(373.15)*1e6 atol=0.0001
@test 0.0640 ≈ Psychro.dBaa(473.15)*1e6 atol=0.0001

# Caaa
@test 2267.0 ≈ Psychro.Caaa(173.15)*1e12 atol=1.0
@test 1409.0 ≈ Psychro.Caaa(273.15)*1e12 atol=1.0
@test 1202.0 ≈ Psychro.Caaa(373.15)*1e12 atol=1.0
@test 1139.0 ≈ Psychro.Caaa(473.15)*1e12 atol=1.0

# dCaaa/dT
@test -18.00 ≈ Psychro.dCaaa(173.15)*1e12 atol=0.01
@test -3.65 ≈ Psychro.dCaaa(273.15)*1e12 atol=0.01
@test -1.06 ≈ Psychro.dCaaa(373.15)*1e12 atol=0.01
@test -0.34 ≈ Psychro.dCaaa(473.15)*1e12 atol=0.01

# Baw
@test -93.3 ≈ Psychro.Baw(173.15)*1e6 atol=0.1
@test -36.4 ≈ Psychro.Baw(273.15)*1e6 atol=0.1
@test -14.5 ≈ Psychro.Baw(373.15)*1e6 atol=0.1
@test -3.1 ≈ Psychro.Baw(473.15)*1e6 atol=0.1

# dBaw/dT
@test 1.01 ≈ Psychro.dBaw(173.15)*1e6 atol=0.01
@test 0.318 ≈ Psychro.dBaw(273.15)*1e6 atol=0.001
@test 0.151 ≈ Psychro.dBaw(373.15)*1e6 atol=0.001
@test 0.0869 ≈ Psychro.dBaw(473.15)*1e6 atol=0.0001

# Caaw
@test 1023 ≈ Psychro.Caaw(173.15)*1e12 atol=1.0
@test 861 ≈ Psychro.Caaw(273.15)*1e12 atol=1.0
@test 696 ≈ Psychro.Caaw(373.15)*1e12 atol=1.0
@test 627 ≈ Psychro.Caaw(473.15)*1e12 atol=1.0

# dCaaw/dT
@test  5.56 ≈ Psychro.dCaaw(173.15)*1e12 atol=0.01
@test -2.44 ≈ Psychro.dCaaw(273.15)*1e12 atol=0.01
@test -1.02 ≈ Psychro.dCaaw(373.15)*1e12 atol=0.01
@test -0.457 ≈ Psychro.dCaaw(473.15)*1e12 atol=0.001

# Caww
@test -2.0e7 ≈ Psychro.Caww(173.15)*1e12 atol=0.1e7
@test -2.2e5 ≈ Psychro.Caww(273.15)*1e12 atol=0.1e5
@test -3.0e4 ≈ Psychro.Caww(373.15)*1e12 atol=0.1e4
@test -8.4e3 ≈ Psychro.Caww(473.15)*1e12 atol=0.1e3

# dCaww/dT
@test 1.61e6 ≈ Psychro.dCaww(173.15)*1e12 atol=0.01e6
@test 6.05e3 ≈ Psychro.dCaww(273.15)*1e12 atol=0.01e3
@test    456 ≈ Psychro.dCaww(373.15)*1e12 atol=1.0
@test   86.9 ≈ Psychro.dCaww(473.15)*1e12 atol=0.1
# Note that on table 1, the last line is 8.69 instead of 86.9. Probably a typo in the paper.

# Specific volume of dry air
# Table 3
P = 0.1e6
@test 0.49510 ≈ Psychro.volumeair(173.15, P) atol=0.00003
@test 0.78364 ≈ Psychro.volumeair(273.15, P) atol=0.0001
@test 1.0713 ≈ Psychro.volumeair(373.15, P) atol=0.0001
@test 1.3586 ≈ Psychro.volumeair(473.15, P) atol=0.0001

P = 0.5e6
@test 0.09747 ≈ Psychro.volumeair(173.15, P) atol=0.00003
@test 0.15637 ≈ Psychro.volumeair(273.15, P) atol=0.00003
@test 0.21436 ≈ Psychro.volumeair(373.15, P) atol=0.00003
@test 0.27207 ≈ Psychro.volumeair(473.15, P) atol=0.00003

P = 1e6
@test 0.04775 ≈ Psychro.volumeair(173.15, P) atol=0.00003
@test 0.07797 ≈ Psychro.volumeair(273.15, P) atol=0.00003
@test 0.10725 ≈ Psychro.volumeair(373.15, P) atol=0.00003
@test 0.13625 ≈ Psychro.volumeair(473.15, P) atol=0.00003

P = 5e6
@test 0.00795 ≈ Psychro.volumeair(173.15, P) atol=0.00003
@test 0.01533 ≈ Psychro.volumeair(273.15, P) atol=0.00003
@test 0.02162 ≈ Psychro.volumeair(373.15, P) atol=0.00003
@test 0.02763 ≈ Psychro.volumeair(473.15, P) atol=0.00003

# Specific enthalpy of dry air
# Table 3
P = 0.1e6
@test -100.63 ≈ Psychro.enthalpyair(173.15, P)/1000 atol=0.01
@test    0.004≈ Psychro.enthalpyair(273.15, P)/1000 atol=0.01
@test  100.79 ≈ Psychro.enthalpyair(373.15, P)/1000 atol=0.01
@test  202.55 ≈ Psychro.enthalpyair(473.15, P)/1000 atol=0.01

P = 0.5e6
@test -103.17 ≈ Psychro.enthalpyair(173.15, P)/1000 atol=0.01
@test   -1.10 ≈ Psychro.enthalpyair(273.15, P)/1000 atol=0.01
@test  100.26 ≈ Psychro.enthalpyair(373.15, P)/1000 atol=0.01
@test  202.30 ≈ Psychro.enthalpyair(473.15, P)/1000 atol=0.01

P = 1e6
@test -106.42 ≈ Psychro.enthalpyair(173.15, P)/1000 atol=0.01
@test   -2.48 ≈ Psychro.enthalpyair(273.15, P)/1000 atol=0.01
@test   99.60 ≈ Psychro.enthalpyair(373.15, P)/1000 atol=0.01
@test  202.00 ≈ Psychro.enthalpyair(473.15, P)/1000 atol=0.01

P = 5e6
@test -135.55 ≈ Psychro.enthalpyair(173.15, P)/1000 atol=0.01
@test  -13.15 ≈ Psychro.enthalpyair(273.15, P)/1000 atol=0.01
@test   94.63 ≈ Psychro.enthalpyair(373.15, P)/1000 atol=0.01
@test  199.8 ≈ Psychro.enthalpyair(473.15, P)/1000 atol=0.1

# Specific entropy of dry air
P = 0.1e6
@test -0.4550 ≈ Psychro.entropyair(173.15, P)/1000 atol=0.0001
@test  0.0038 ≈ Psychro.entropyair(273.15, P)/1000 atol=0.0001
@test  0.3182 ≈ Psychro.entropyair(373.15, P)/1000 atol=0.0001
@test  0.5597 ≈ Psychro.entropyair(473.15, P)/1000 atol=0.0001

# Table 4, appendix of reference [2]
# Testing the volume of saturated moist air
volmoist(Tk, P) = Psychro.volumemoist(Tk, P, Psychro.efactor(Tk,P)*Psychro.Pws(Tk)/P)
P = 0.1e6
@test 0.49510 ≈ volmoist(173.15, P) atol=0.00005
@test 0.78846 ≈ volmoist(273.15, P) atol=0.00005
@test 3.4974 ≈ volmoist(363.15, P) atol=0.0002
P = 0.5e6
@test 0.097466 ≈ volmoist(173.15, P) atol=0.000006
@test 0.15657 ≈ volmoist(273.15, P) atol=0.00005
@test 0.24273 ≈ volmoist(363.15, P) atol=0.00003
@test 0.26904 ≈ volmoist(373.15, P) atol=0.00003
P=1e6
@test 0.047752 ≈ volmoist(173.15, P) atol=0.000006
@test 0.078022 ≈ volmoist(273.15, P) atol=0.000005
@test 0.11226 ≈ volmoist(363.15, P) atol=0.00001
@test 0.11942 ≈ volmoist(373.15, P) atol=0.00001
@test 0.23281 ≈ volmoist(423.15, P) atol=0.00002
@test 0.32711 ≈ volmoist(433.15, P) atol=0.00002
@test 0.61489 ≈ volmoist(443.15, P) atol=0.00005

P=5e6
@test 0.007951 ≈ volmoist(173.15, P) atol=0.000006
@test 0.015331 ≈ volmoist(273.15, P) atol=0.000005
@test 0.021310 ≈ volmoist(363.15, P) atol=0.00001
@test 0.022073 ≈ volmoist(373.15, P) atol=0.00001
@test 0.02727  ≈ volmoist(423.15, P) atol=0.00001
@test 0.02885  ≈ volmoist(433.15, P) atol=0.00001
@test 0.03077  ≈ volmoist(443.15, P) atol=0.00001
@test 0.03315  ≈ volmoist(453.15, P) atol=0.00001
@test 0.0402  ≈ volmoist(473.15, P) atol=0.0001

hmoist(Tk, P) = Psychro.enthalpymoist(Tk, P, Psychro.efactor(Tk,P)*Psychro.Pws(Tk)/P)/1000
P = 0.1e6
@test -100.627 ≈ hmoist(173.15, P) atol=0.005
@test    9.602 ≈ hmoist(273.15, P) atol=0.005
@test   4035.0 ≈ hmoist(363.15, P) atol=0.3 # In the paper the number is 4305!

P = 0.5e6
@test -103.168 ≈ hmoist(173.15, P) atol=0.005
@test    0.831 ≈ hmoist(273.15, P) atol=0.005
@test   365.55 ≈ hmoist(363.15, P) atol=0.03
@test   533.27 ≈ hmoist(373.15, P) atol=0.03
@test   35957 ≈ hmoist(423.15, P) atol=2.0

P = 1e6
@test -106.415 ≈ hmoist(173.15, P) atol=0.005
@test   -1.498 ≈ hmoist(273.15, P) atol=0.005
@test   217.77 ≈ hmoist(363.15, P) atol=0.03
@test   293.17 ≈ hmoist(373.15, P) atol=0.03
@test   1788.8 ≈ hmoist(423.15, P) atol=0.5
@test   3113.3 ≈ hmoist(433.15, P) atol=0.5
@test   7206.0 ≈ hmoist(443.15, P) atol=2.0

P = 5e6
@test -135.547 ≈ hmoist(173.15, P) atol=0.005
@test  -12.929 ≈ hmoist(273.15, P) atol=0.005
@test   109.89 ≈ hmoist(363.15, P) atol=0.03
@test   132.42 ≈ hmoist(373.15, P) atol=0.03
@test   347.3 ≈ hmoist(423.15, P) atol=0.5
@test   428.8 ≈ hmoist(433.15, P) atol=0.5
@test   534.1 ≈ hmoist(443.15, P) atol=0.5
@test   672.0 ≈ hmoist(453.15, P) atol=1.0
@test   1113.0 ≈ hmoist(473.15, P) atol=1.0

# Entropy of saturated moist air
# Table 4 in the appendix of reference [2].
smoist(Tk, P) = Psychro.entropymoist(Tk, P, Psychro.efactor(Tk,P)*Psychro.Pws(Tk)/P)/1000
P = 0.1e6
@test  -0.45499 ≈ smoist(173.15, P) atol=0.00003
@test   0.04070 ≈ smoist(273.15, P) atol=0.00003
@test  11.7287 ≈ smoist(363.15, P) atol=0.0005

P = 0.5e6
@test  -0.92718 ≈ smoist(173.15, P) atol=0.00005
@test  -0.45417 ≈ smoist(273.15, P) atol=0.00005
@test   0.64517 ≈ smoist(363.15, P) atol=0.0005
@test   1.1102 ≈ smoist(373.15, P) atol=0.0005
@test  90.01 ≈ smoist(423.15, P) atol=0.05

P = 1e6
@test  -1.13929 ≈ smoist(173.15, P) atol=0.00005 # In the table the signal is +
@test  -0.66104 ≈ smoist(273.15, P) atol=0.00005
@test   0.00805 ≈ smoist(363.15, P) atol=0.00005
@test   0.21676 ≈ smoist(373.15, P) atol=0.0005
@test   4.0478  ≈ smoist(423.15, P) atol=0.0005
@test   7.3103  ≈ smoist(433.15, P) atol=0.0005
@test   17.221  ≈ smoist(443.15, P) atol=0.005

P = 5e6
@test  -1.72380 ≈ smoist(173.15, P) atol=0.0001 # In the table the signal is +
@test  -1.15924 ≈ smoist(273.15, P) atol=0.00005
@test  -0.77519 ≈ smoist(363.15, P) atol=0.00005
@test  -0.71331 ≈ smoist(373.15, P) atol=0.00005
@test  -0.1629  ≈ smoist(423.15, P) atol=0.0002
@test   0.0364  ≈ smoist(433.15, P) atol=0.0002
@test   0.2901  ≈ smoist(443.15, P) atol=0.0002
@test   0.6185  ≈ smoist(453.15, P) atol=0.0002
@test   1.642  ≈ smoist(473.15, P) atol=0.002


text:
"""
supported 'dataType' attributes for dataSets
"CT" - will lead to classical display of CT image
"boolLabel" - binary color display (either some color or none)
"multiDiscreteLabel" - diffrent discrete colors assigned to multiple discrete labels
"contLabel" - continous color  display for example float valued probabilistic label
"PET" - continous color  display for example float valued
"manualModif" - manually modified array - using annotation functionality of a viewer

"""

module visualizationFromHdf5

export getGroupOrCreate,getArrByName,loadFromHdf5Prim,getSomeColor,openHDF5,calculateAndDisplay,writeGroupAttribute

import MedEye3d
import MedEye3d.ForDisplayStructs
import MedEye3d.ForDisplayStructs.TextureSpec
using ColorTypes
import MedEye3d.SegmentationDisplay
using HDF5
import MedEye3d.DataStructs.ThreeDimRawDat
import MedEye3d.DataStructs.DataToScrollDims
import MedEye3d.DataStructs.FullScrollableDat
import MedEye3d.ForDisplayStructs.KeyboardStruct
import MedEye3d.ForDisplayStructs.MouseStruct
import MedEye3d.ForDisplayStructs.ActorWithOpenGlObjects
import MedEye3d.OpenGLDisplayUtils
import MedEye3d.DisplayWords.textLinesFromStrings
import MedEye3d.StructsManag.getThreeDims

using ..distinctColorsSaved

using MedEval3D
using MedEval3D.BasicStructs
using MedEval3D.MainAbstractions

"""
simply open dataset
"""
function openHDF5(pathToHDF5)
  return  h5open(pathToHDF5, "r+")

end

"""
Writes group attribute - simple wrapper
group - object representing the group
attrName - name of the attribute
value - value of the attribute
"""
function writeGroupAttribute(fid,groupName, attrName, value)
  group = fid[groupName]
  write_attribute(group, attrName, value)
end

"""
return the HDF5 group if it is not already present in dataset it creates it
  fid - HDF5 object
  groupName - string representing group of intrest
"""
function getGroupOrCreate(fid, groupName)
  if(!haskey(fid, groupName))
    return create_group(fid, groupName)
  end
  #in case it is already created
  return fid[groupName]

end

"""
get some color from listOfColors or if those are already used some random color from
listOfColorUsed- boolean array marking which colors were already used from listOfColors
"""
function getSomeColor(listOfColorUsed)
  if(sum(listOfColorUsed)<18)
  for i in 1:18
    if(!listOfColorUsed[i])
      listOfColorUsed[i]=true;
      tupl= listOfColors[i]
      return RGB(tupl[1]/255,tupl[2]/255,tupl[3]/255)
    end#if
  end#for
else
   #if we are here it means that no more colors from listOfColors is available - so we need to take some random color from bigger list
   tupl = longColorList[rand(Int,1:255)]

   return  RGB(tupl[1]/255,tupl[2]/255,tupl[3]/255)
end
end#getSomeColor

"""
becouse Julia arrays is column wise contiguus in memory and open GL expects row wise we need to rotate and flip images
pixels - 3 dimensional array of pixel data
"""
function permuteAndReverse(pixels)

    pixels=  permutedims(pixels, (3,2,1))
    sizz=size(pixels)
    for i in 1:sizz[2]
        for j in 1:sizz[3]
            pixels[:,i,j] =  reverse(pixels[:,i,j])
        end#
    end#
    return pixels
  end#permuteAndReverse

function onlyPermute(pixels)
    pixels=  permutedims(pixels, (3,2,1))

    sizz=size(pixels)
    for i in 1:sizz[1]
        for j in 1:sizz[3]
            pixels[i,:,j] =  reverse(pixels[i,:,j])
        end#
    end#

    return pixels
  end#permuteAndReverse

"""
loading image from HDF5 dataset
fid- objct managing HDF5 file
patienGroupName- tell us  the string that is a name of the HDF5 group with all data needed for a given patient
fractionOfMainIm - effectively will set how much space is left for text display
addTextSpecs - additional texture specifications required to define and save some of the masks
"""
function loadFromHdf5Prim(fid, patienGroupName::String
  , addTextSpecs::Vector{TextureSpec},listOfColorUsed,fractionOfMainIm::Float32= Float32(0.8) )

#marks what colors are already used
  group = fid[patienGroupName]
 #strings holding the arrays holding data about given patient
 imagesMasks = keys(group)

#adding one spot to be able to get manually modifiable mask
 imageSize::Tuple{Int64,Int64,Int64}= (0,0,0)
toAddAddedTextures = 0
addTextureNames::Vector{String}= Vector(undef,length(addTextSpecs))
#making sure that each array will be included
index = 0;
for tex in addTextSpecs
  index+=1
  addTextureNames[index]=tex.name
  if(!haskey(group,tex.name))
    toAddAddedTextures+=1
  end#if
end#for

#preinitialized placeholders
 textureSpecifications::Vector{TextureSpec}=Vector(undef,length(imagesMasks)+toAddAddedTextures)
 tupleVect=Vector(undef,length(imagesMasks)+toAddAddedTextures)

 index = 0;
  for maskName in imagesMasks
    index+=1
    dset = group[maskName]
    dataTypeStr= attributes(dset)["dataType"][]
    #additional arrays are saved and loaded as is
    voxels =  dset[:,:,:]

    # @info " mask name   " maskName
    # @info "voxel dims   " size(voxels)

    typp = eltype(voxels)
    min = attributes(dset)["min"][]
    max = attributes(dset)["max"][]

#initializing texture
    textureSpec =  if(maskName in addTextureNames )
      textureSpec = filter(tex->tex.name==maskName,addTextSpecs  )[1]
    else getDefaultTextureSpec(dataTypeStr,maskName ,index,listOfColorUsed, typp, min, max)
  end
  textureSpec.numb= index
    imageSize= size(voxels)

    textureSpecifications[index]=textureSpec
    #append!( textureSpecifications, textureSpec )
    tupleVect[index] =(maskName,voxels)
  end #for
#and additionally manually modifiable ...

for tex in addTextSpecs
  if( !haskey(group, tex.name) )
    index+=1

    #@info "nnnnnn no  manualModif key "
    textureSpec = tex
    textureSpec.numb= index
    textureSpecifications[index]=textureSpec
    tupleVect[index] =(tex.name,zeros(eltype(tex.minAndMaxValue) ,imageSize))
  end#if
end#for

spacingList = attributes(group)["spacing"][]
spacing=(Int64(spacingList[1]),Int64(spacingList[2]),Int64(spacingList[3]))

datToScrollDimsB= MedEye3d.ForDisplayStructs.DataToScrollDims(
    imageSize=  imageSize
    ,voxelSize=spacing
    ,dimensionToScroll = 3 );

slicesDat= getThreeDims(tupleVect )

#mainLines= textLinesFromStrings([""]);
#supplLines=[];

mainLines= textLinesFromStrings([" "]);
supplLines=map(x->  textLinesFromStrings(["slice $(x)"]), 1:imageSize[3] );

mainScrollDat = FullScrollableDat(dataToScrollDims =datToScrollDimsB
                                 ,dimensionToScroll=1 # what is the dimension of plane we will look into at the beginning for example transverse, coronal ...
                                 ,dataToScroll= slicesDat
                                 ,mainTextToDisp= mainLines
                                 ,sliceTextToDisp=supplLines );

                                 SegmentationDisplay.coordinateDisplay(textureSpecifications
                                 ,fractionOfMainIm ,datToScrollDimsB ,1000);

                                 SegmentationDisplay.passDataForScrolling(mainScrollDat);

return mainScrollDat

 end #loadFromHdf5Prim

 """
 given datatype string will return appropriate default configuration for image display
 """
function getDefaultTextureSpec(dataTypeStr::String,maskName::String ,index::Int,listOfColorUsed
      , typp, min, max)::TextureSpec

  if(dataTypeStr=="CT")
    return TextureSpec{typp}(
      name= maskName,
      isMainImage = true,
      minAndMaxValue= typp.([0,100]))

  elseif(dataTypeStr=="boolLabel")
    return TextureSpec{typp}(
      name = maskName,
      color =getSomeColor(listOfColorUsed)
      ,minAndMaxValue= typp.([min,max])
     )

  elseif(dataTypeStr=="multiDiscreteLabel")

  elseif(dataTypeStr=="contLabel")

  elseif(dataTypeStr=="PET")

  end

end

function getArrByName(name::String,mainScrollDat )

return  filter((it)->it.name==name ,mainScrollDat.dataToScroll)[1].dat

end #getArrByName

function calculateAndDisplay(preparedDict,mainScrollDat::FullScrollableDat
  , conf::ConfigurtationStruct,numberToLookFor,goldArr,algoOutputGPU)
  res= calcMetricGlobal(preparedDict,conf,goldArr,algoOutputGPU,numberToLookFor)
  append!(mainScrollDat.mainTextToDisp, textLinesFromStrings(giveStringsFromResultMetrics(res,conf)) )

end

"""
supplied ResultMetrics struct will return list of strings with results
"""
function giveStringsFromResultMetrics(res,conf::ConfigurtationStruct)::Vector{String}
  output = []
  if(conf.dice)
    append!(output,["dice $(res.dice)"])
  end
  if(conf.jaccard)
    append!(output,["jaccard $(res.jaccard)"])
  end
  if(conf.gce)
    append!(output,["gce $(res.gce)"])
  end
  if(conf.vol)
    append!(output,["vol $(res.vol)"])
  end
  if(conf.randInd)
    append!(output,["randInd $(res.randInd)"])
  end
  if(conf.ic)
    append!(output,["ic $(res.ic)"])
  end
  if(conf.kc)
    append!(output,["kc $(res.kc)"])
  end
  if(conf.mi)
    append!(output,["mi $(res.mi)"])
  end
  if(conf.vi)
    append!(output,["vi $(res.vi)"])
  end
  if(conf.md)
    append!(output,["md $(res.md)"])
  end
  if(conf.hd)
    append!(output,["hd $(res.hd)"])
  end

  return output

end

end#end


text:
using JuMP

x0 = rand(9)

function foo(x...)
    sum(collect(x).^2)
end

function jump_energy(x0)
    model = Model()
    @defVar(model, x[i=1:length(x0)], start=x0[i])
    registerNLFunction(:foo, length(x0), foo, autodiff=true)
    @setNLObjective(model, Min, foo(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9]))
    solve(model)
    @show getValue(x)
end

jump_energy(x0)


text:
using BioBackgroundModels, Serialization

hmm_output = "/bench/PhD/NGS_binaries/BBM/survey_chains"
sample_output = "/bench/PhD/NGS_binaries/BBM/survey_samples"

survey_folders = "/bench/PhD/NGS_binaries/BBM/survey_folders"

chains=deserialize(hmm_output)
sample_dfs = deserialize(sample_output)
training_sets, test_sets = split_obs_sets(sample_dfs)

report_folders=generate_reports(chains, test_sets)
serialize(survey_folders, report_folders)


text:
# currently supported packages
const SUPPORTED_EXPONENTIAL_PACKAGES = [:ExponentialUtilities, :Expokit]

# =============
# Global option
# =============

global exponential_backend = missing  # global state of the exponential backend

function set_exponential_backend!(backend::Module)
    global exponential_backend = Val(Symbol(backend))
end

function get_exponential_backend()
    ismissing(exponential_backend) && error("no exponential backend is loaded; " *
        "load one of these packages: $SUPPORTED_EXPONENTIAL_PACKAGES")
    return exponential_backend
end

# ====================
# ExponentialUtilities
# ====================

function load_exponentialutilities()
return quote

if ismissing(exponential_backend)
    set_exponential_backend!(ExponentialUtilities)
end

end end  # quote / load_exponentialutilities

function _expmv(::Val{:ExponentialUtilities}, t, A, b)
    return ExponentialUtilities.expv(t, A, b)
end

# =======
# Expokit
# =======

function load_expokit()
return quote

if ismissing(exponential_backend)
    set_exponential_backend!(Expokit)
end

end end  # quote / load_expokit

function _expmv(::Val{:Expokit}, t, A, b::Vector)
    return Expokit.expmv(t, A, b)
end

function _expmv(backend::Val{:Expokit}, t, A, b)
    b = Vector(b)  # Expokit requires a dense vector
    return _expmv(backend, t, A, b)
end


text:
cd(@__DIR__); import Pkg; Pkg.activate(".")
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
Pkg.develop(path="$(homedir())/github/SimpleSDMLayers.jl")
using Revise
using SimpleSDMLayers
using Test

## Step 6: Contiguous hcat/vcat ####
using Test

# Base layer
l1 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=0.0, right=10.0, bottom=0.0, top=10.0)
# Non contiguous layers
l2 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=0.0, right=10.0, bottom=20.0, top=30.0)
l3 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=20.0, right=30.0, bottom=0.0, top=10.0)

vcat(l2, l1)
hcat(l1, l3)

l4 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=0.0, right=10.0, bottom=10.0, top=20.0)
vcat(l4, l1)
vcat(l1, l4)

## Step 7: No replace! in geotiff ####
import SimpleSDMLayers._prepare_layer_for_burnin
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
layer = SimpleSDMPredictor(WorldClim, BioClim, 1)

array = layer.grid
replace!(array, nothing => NaN)
layer.grid # NaN

array = replace(layer.grid, nothing => NaN)
layer.grid # NaN


text:
#==============================================================

    MAKx.jl

    different formula for each cycle (dfc)

    write functions to fit MAKx (MAK2 and MAK3) model here,
    which will be called in fit_amplification_model()
    in amplification.jl

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

import DataStructures.OrderedDict;
import JuMP: Model, @variable, @constraint, @NLconstraint, @NLobjective,
    solve, getvalue, getobjectivevalue

## start values
const MAK_d0_START  = 0  ## 0 good, 1 bad
const k_START       = 10 ## used: 10 better, 2 good, 1e-10 bad, 1 bad

## function definitions

pred_from_d_nm1(::Union{Type{Val{MAK2}}, Type{Val{MAK3}}}, d_nm1::Real, k::Real) =
    d_nm1 + k * log(1 + d_nm1 / k)

function pred_from_cycs( ## 0.7to1.2e-5 sec for 40 cycles on PC
    ::Type{Val{MAK2}},
    cycs    ::AbstractVector,
    fb      ::Real,
    d0      ::Real,
    k       ::Real
)
    max_cyc = maximum(cycs)
    pred_ds = [Float_T(d0)]
    i = 1
    while i <= max_cyc
        push!(pred_ds, pred_from_d_nm1(MAK2(), pred_ds[i], k))
        i += 1
    end
    return fb + pred_ds[2:end][map(Int_T, cycs)]
end ## pred_from_cycs(::Type{Val{MAK2}})

function pred_from_cycs( ## sec for 40 cycles on PC
    ::Type{Val{MAK3}},
    cycs    ::AbstractVector,
    fb      ::Real,
    bl_k    ::Real,
    d0      ::Real,
    k       ::Real
)
    max_cyc = maximum(cycs)
    pred_ds = [Float_T(d0)]
    i = 1
    while i <= max_cyc
        push!(pred_ds, pred_from_d_nm1(MAK2(), pred_ds[i], k))
        i += 1
    end
    return fb + bl_k * cycs .+ pred_ds[2:end][map(Int_T, cycs)]
end ## pred_from_cycs(::Type{Val{MAK3}})

function fit(
    ::Type{Val{MAK2}},
    cycs        ::AbstractVector, ## continous integers or not
    obs_fluos   ::AbstractVector,
    wts         ::AbstractVector =ones(length(obs_fluos));
    kwargs_Model... # argument for `JuMP.Model`
)
    ## find approximate `max_d` by `finite_diff` from "shared.jl"
    d_vec = finite_diff(cycs, obs_fluos; nu=1) ## should use `nu=2` per Boggy 2010 paper (because MAK2 is no longer valid at max_d1), but use `nu=1` to increase the number of data points used for fit and seems to provide better fitting results
    max_d, max_d_idx = findmax(d_vec)
    idc2fit = 1:max_d_idx
    obs2fit = obs_fluos[idc2fit]
    cycs2fit = 1:cycs[max_d_idx]

    jmp_model = Model(;kwargs_Model...)

    @variable(jmp_model, fb, start=minimum(obs2fit))
    @variable(jmp_model, d0 >= 0, start=MAK_d0_START)
    @variable(jmp_model, k >= 1e-10, start=k_START)
    # @variable(jmp_model, bl_k, start=0)
    @variable(jmp_model, f[cycs2fit])
    @variable(jmp_model, d[cycs2fit]) ## change_a1
    # @variable(jmp_model, d[cycs2fit] >= 0) ## change_a2 # didn't work as well as "change_a1", see "20170312_0302_ip137_exp187_ch1_mak2_ylims_obs_*.png"

    @constraint(jmp_model, f_constr[cyc in cycs2fit], f[cyc] == fb + d[cyc])
    @NLconstraint(jmp_model, d_constr_01, d[1] == d0 + k * log(1 + d0 / k))
    @NLconstraint(
        jmp_model,
        d_constr_2p[cyc in cycs2fit[2:end]],
        d[cyc] == d[cyc-1] + k * log(1 + d[cyc-1] / k)
    )

    @NLobjective(
        jmp_model, Min,
        sum((f[cycs[idx]] - obs_fluos[idx]) ^ 2 for idx in idc2fit)
    )

    # @NLobjective(
    #     jmp_model, Min,
    #     @eval $(parse(get_mak2_rsq_str(obs_fluos[1:max_d_idx]))) ## `OutOfMemoryError()`
    # )

    status = solve(jmp_model)
    coef_strs = ["fb", "d0", "k"]
    coefs = map(getvalue, [fb, d0, k])
    obj_val = getobjectivevalue(jmp_model)

    return MAK2Fit(
        max_d_idx,
        coef_strs,
        coefs,
        status,
        obj_val,
        jmp_model,
        # init_coefs
    )
end ## fit(::Type{Val{MAK2}})

function fit(
    ::Type{Val{MAK3}},
    cycs        ::AbstractVector, ## continous integers or not
    obs_fluos   ::AbstractVector,
    wts         ::AbstractVector =ones(length(obs_fluos));
    kwargs_Model... ## argument for `JuMP.Model`
)
    ## find approximate `max_d` by `finite_diff` from "shared.jl"
    d_vec = finite_diff(cycs, obs_fluos; nu=1) ## should use `nu=2` per Boggy 2010 paper (because MAK2 is no longer valid at max_d1), but use `nu=1` to increase the number of data points used for fit and seems to provide better fitting results
    max_d, max_d_idx = findmax(d_vec)
    idc2fit = 1:max_d_idx
    cycs2fit = 1:cycs[max_d_idx]

    ## fit a linear model to the estimated baseline portion of the curve
    d2_vec = finite_diff(cycs, obs_fluos; nu=2)
    max_d2, max_d2_idx = findmax(d2_vec)
    idc2fit_4bl = 1:(max(1, max_d2_idx - 1))
    fb_start, bl_k_start = linreg(cycs[idc2fit_4bl], obs_fluos[idc2fit_4bl])

    jmp_model = Model(;kwargs_Model...)

    @variable(jmp_model, fb, start=fb_start)
    @variable(jmp_model, bl_k, start=bl_k_start)
    @variable(jmp_model, d0 >= 0, start=MAK_d0_START)
    @variable(jmp_model, k >= 1e-10, start=k_START)
    # @variable(jmp_model, bl_k, start=0)
    @variable(jmp_model, f[cycs2fit])
    @variable(jmp_model, d[cycs2fit])

    @constraint(jmp_model,
        f_constr[cyc in cycs2fit],
        f[cyc] == fb + bl_k * cyc + d[cyc]
    )
    @NLconstraint(jmp_model, d_constr_01, d[1] == d0 + k * log(1 + d0 / k))
    @NLconstraint(
        jmp_model,
        d_constr_2p[cyc in cycs2fit[2:end]],
        d[cyc] == d[cyc-1] + k * log(1 + d[cyc-1] / k)
    )

    @NLobjective(
        jmp_model, Min,
        sum((f[cycs[idx]] - obs_fluos[idx]) ^ 2 for idx in idc2fit)
    )

    # @NLobjective(
    #     jmp_model, Min,
    #     @eval $(parse(get_mak2_rsq_str(obs_fluos[1:max_d_idx]))) ## `OutOfMemoryError()`
    # )

    status = solve(jmp_model)
    coef_strs = ["fb", "bl_k", "d0", "k"]
    coefs = map(getvalue, [fb, bl_k, d0, k])
    obj_val = getobjectivevalue(jmp_model)

    return MAK3Fit(
        max_d_idx,
        fb_start,
        bl_k_start,
        coef_strs,
        coefs,
        status,
        obj_val,
        jmp_model,
        # init_coefs
    )
end ## fit(::Type{Val{MAK3}})

## need `@eval`, not used
#
#
# const d_nm1_STR = "d_nm1"
# const k_STR = "k"
# const d0_STR = "d0"
# const fb_STR = "fb"
#
#
## associated with "get_mak2_pff"
# function get_mak2_psf(
#     num_cycs::Integer;
#     d_nm1_str::String=d_nm1_STR,
#     k_str::String=k_STR,
#     d0_str::String=d0_STR
#     )
#     pred_str_d_nm1_n = "($d_nm1_str + $k_str * log(1 + $d_nm1_str / $k_str))"
#     pred_strs_d_w0 = [d0_str]
#     wo0_range = 2:num_cycs+1
#     for cyc_n in wo0_range
#         push!(
#             pred_strs_d_w0,
#             replace(pred_str_d_nm1_n, d_nm1_str, pred_strs_d_w0[cyc_n - 1])
#         )
#     end # for
#     pred_strs_fluo = map(pred_strs_d_w0[wo0_range]) do pred_str_d
#         "fb + $pred_str_d"
#     end
#     return pred_strs_fluo
# end ## get_mak2_psf
#
#
## `OutOfMemoryError()` upon macro expansion
# function get_mak2_pff( # return an array whose each element is a prediction function for each cycle.
#     num_cycs::Integer;
#     d_nm1_str::String=d_nm1_STR,
#     k_str::String=k_STR,
#     d0_str::String=d0_STR,
#     fb_str::String=fb_STR
#     )
#     pred_strs_fluo = get_mak2_psf(num_cycs;
#         d_nm1_str=d_nm1_str,
#         k_str=k_str,
#         d0_str=d0_str
#     )
#     return map(1:num_cycs) do cyc_n
#         func_expr = parse("function pred_$cyc_n($fb_str::Real, $d0_str::Real, $k_str::Real) return $(pred_strs_fluo[cyc_n]) end")
#         @eval $func_expr # `OutOfMemoryError()`
#     end ## do pred_str_fluo
# end
#
#
## for `@NLobjective(jmp_model, Min, @eval ...)` from `fit_mak2`
# function get_mak2_rsq_str(
#     obs_fluos::AbstractVector;
#     kwargs_gmp::Associative=OrderedDict(), # keyword arguments passed onto `get_mak2_psf`
#     )
#     num_cycs = length(obs_fluos)
#     pred_strs_fluo = get_mak2_psf(num_cycs; kwargs_gmp...)
#     rsq_str = join(map(1:num_cycs) do cyc_n
#         "($(pred_strs_fluo[cyc_n]) - $(obs_fluos[cyc_n])) ^ 2"
#     end, " + ") ## do cyc_n
#     return rsq_str
# end ## get_mak2_rsq_str

## for later use
# pred_funcs = get_mak2_pff(num_cycs)
# preds = map(func -> func(fb, d0, k), funcs)

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
export MonteCarloLearner,
       FIRST_VISIT,
       EVERY_VISIT,
       NO_SAMPLING,
       ORDINARY_IMPORTANCE_SAMPLING,
       WEIGHTED_IMPORTANCE_SAMPLING

using StatsBase: countmap

abstract type AbstractVisitType end
struct FirstVisit <: AbstractVisitType end
struct EveryVisit <: AbstractVisitType end
const FIRST_VISIT = FirstVisit()
const EVERY_VISIT = EveryVisit()

abstract type AbstractSamplingStyle end
struct NoSampling <: AbstractSamplingStyle end
struct OrdinaryImportanceSampling <: AbstractSamplingStyle end
struct WeightedImportanceSampling <: AbstractSamplingStyle end
const NO_SAMPLING = NoSampling()
const ORDINARY_IMPORTANCE_SAMPLING = OrdinaryImportanceSampling()
const WEIGHTED_IMPORTANCE_SAMPLING = WeightedImportanceSampling()

"""
    MonteCarloLearner(;kwargs...)

# Fields

- `approximator`::[`AbstractApproximator`](@ref),
"""
struct MonteCarloLearner{T,A<:AbstractApproximator,R,S} <: AbstractLearner
    approximator::A
    γ::Float64
    α::Float64
    returns::R

    MonteCarloLearner(
        ;
        approximator::A,
        γ = 1.0,
        α = 1.0,
        kind = FIRST_VISIT,
        sampling = NO_SAMPLING,
        returns = CachedSampleAvg{Float64}(),
    ) where {A} =
        new{typeof(kind),A,typeof(returns),typeof(sampling)}(approximator, γ, α, returns)
end

VisitStyle(::MonteCarloLearner{T}) where T = T
SamplingStyle(::MonteCarloLearner{T,A,R,S}) where {T,A,R,S} = S

(learner::MonteCarloLearner)(obs) = learner.approximator(get_state(obs))
(learner::MonteCarloLearner)(obs, a) = learner.approximator(get_state(obs))[a]

RLBase.update!(learner::MonteCarloLearner, experience::NamedTuple) = update!(learner, VisitStyle(learner), ApproximatorStyle(learner.approximator), SamplingStyle(learner), experience)

function extract_experience(t::AbstractTrajectory, learner::MonteCarloLearner)
    # only extract & update at the end of an episode
    if length(t) > 0 && get_trace(t, :terminal)[end]
        (
            states=get_trace(t, :state),
            actions = get_trace(t, :action),
            rewards = get_trace(t, :reward),
        )
    else
        nothing
    end
end

function RLBase.update!(
    learner::MonteCarloLearner,
    ::Type{FirstVisit},
    ::VApproximator,
    ::Type{NoSampling},
    transitions::NamedTuple{(:states, :actions, :rewards)},
)
    states, rewards = transitions.states, transitions.rewards
    V, γ, α, Returns, G, T = learner.approximator,
        learner.γ,
        learner.α,
        learner.returns,
        0.0,
        length(states)
    seen_states = countmap(states)

    for t = T:-1:1
        S, R = states[t], rewards[t]
        G = γ * G + R
        if seen_states[S] == 1  # first visit
            update!(V, S => α * (Returns(S, G) - V(S)))
            delete!(seen_states, S)
        else
            seen_states[S] -= 1
        end
    end
end

function RLBase.update!(
    learner::MonteCarloLearner,
    ::Type{FirstVisit},
    ::VApproximator,
    ::Type{OrdinaryImportanceSampling},
    transitions::NamedTuple{(:states, :actions, :rewards, :weights)},
)
    states, rewards, weights = transitions.states, transitions.rewards, transitions.weights
    V, γ, α, Returns, G, ρ, T = learner.approximator,
        learner.γ,
        learner.α,
        learner.returns,
        0.0,
        1.0,
        length(states)
    seen_states = countmap(states)

    for t = T:-1:1
        S, R = states[t], rewards[t]
        G = γ * G + R
        ρ *= weights[t]
        if seen_states[S] == 1  # first visit
            update!(V, S => α * (Returns(S, ρ * G) - V(S)))
            delete!(seen_states, S)
        else
            seen_states[S] -= 1
        end
    end
end

function RLBase.update!(
    learner::MonteCarloLearner,
    ::Type{FirstVisit},
    ::VApproximator,
    ::Type{WeightedImportanceSampling},
    transitions::NamedTuple{(:states, :actions, :rewards, :weights)},
)
    states, rewards, weights = transitions.states, transitions.rewards, transitions.weights
    V, γ, α, (G_cached, ρ_cached), G, ρ, T = learner.approximator,
        learner.γ,
        learner.α,
        learner.returns,
        0.0,
        1.0,
        length(states)
    seen_states = countmap(states)

    for t = T:-1:1
        S, R = states[t], rewards[t]
        G = γ * G + R
        ρ *= weights[t]
        if seen_states[S] == 1  # first visit
            numerator = G_cached(S, ρ * G)
            denominator = ρ_cached(S, ρ)
            val = denominator == 0 ? 0 : numerator / denominator
            update!(V, S => α * (val - V(S)))
            delete!(seen_states, S)
        else
            seen_states[S] -= 1
        end
    end
end

function RLBase.update!(
    learner::MonteCarloLearner,
    ::Type{EveryVisit},
    ::VApproximator,
    ::Type{NoSampling},
    transitions::NamedTuple{(:states, :actions, :rewards)},
)
    states, rewards = transitions.states, transitions.rewards
    α, γ, V, Returns, G = learner.α, learner.γ, learner.approximator, learner.returns, 0.0
    for (s, r) in Iterators.reverse(zip(states, rewards))
        G = γ * G + r
        update!(V, s => α * (Returns(s, G) - V(s)))
    end
end

function RLBase.update!(
    learner::MonteCarloLearner,
    ::Type{FirstVisit},
    ::QApproximator,
    ::Type{NoSampling},
    transitions::NamedTuple{(:states, :actions, :rewards)},
)
    states, actions, rewards = transitions
    α, γ, Q, Returns, G, T = learner.α,
        learner.γ,
        learner.approximator,
        learner.returns,
        0.0,
        length(states)
    seen_pairs = countmap(zip(states, actions))

    for t = T:-1:1
        S, A, R = states[t], actions[t], rewards[t]
        pair = (S, A)
        G = γ * G + R
        if seen_pairs[pair] == 1  # first visit
            update!(Q, pair => α * (Returns(pair, G) - Q(S, A)))
            delete!(seen_pairs, pair)
        else
            seen_pairs[pair] -= 1
        end
    end
end

function RLBase.update!(
    learner::MonteCarloLearner,
    ::Type{EveryVisit},
    ::QApproximator,
    ::Type{NoSampling},
    transitions::NamedTuple{(:states, :actions, :rewards)},
)
    states, actions, rewards = transitions
    α, γ, Q, Returns, G = learner.α, learner.γ, learner.approximator, learner.returns, 0.0
    for (s, a, r) in Iterators.reverse(zip(states, actions, rewards))
        G = γ * G + r
        update!(Q, (s, a) => α * (Returns((s, a), G) - Q(s, a)))
    end
end


text:
using Bridge, Distributions
using Test
#import Bridge: b, σ, a, transitionprob
qu(x) = x*x'
S = SVector{2,Float64}
B = SMatrix{2,2,Float64}([-1 0.1; -0.2 -1])
mu = S([0.2, 0.3])
sigma = SMatrix{2,2,Float64}([-0.212887  0.0687025
  0.193157  0.388997 ])
a = qu(sigma)
@test supnorm(qu(cholupper(a)') - a) < eps()

@test norm(cov(Bridge.mat(S[ cholupper(a)'*randn(S) for i in 1:100])') - Matrix(a)) < 0.1

x = [0.289639+0.649503im, 0.724532+0.963363im]

@test dot(x,x) == Bridge.inner(x)
@test dot(x, 1 .- x) == Bridge.inner(x, 1 .- x)
@test x*x' == Bridge.outer(x)
@test x*(1 .- x)' == Bridge.outer(x, 1 .- x)

A = reshape(collect(1:5*4*3), 5,4,3)
@test A[.., 1] == reshape(1:5*4, 5,4)
A[.., 1] = A[.., 3]
@test A[.., 1] == reshape((3-1) * 5 * 4 .+ (1:5*4), 5, 4)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""ConfigMap holds configuration data for pods to consume.

    IoK8sApiCoreV1ConfigMap(;
        apiVersion=nothing,
        binaryData=nothing,
        data=nothing,
        kind=nothing,
        metadata=nothing,
    )

    - apiVersion::String : APIVersion defines the versioned schema of this representation of an object. Servers should convert recognized schemas to the latest internal value, and may reject unrecognized values. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#resources
    - binaryData::Dict{String, Vector{UInt8}} : BinaryData contains the binary data. Each key must consist of alphanumeric characters, '-', '_' or '.'. BinaryData can contain byte sequences that are not in the UTF-8 range. The keys stored in BinaryData must not overlap with the ones in the Data field, this is enforced during validation process. Using this field will require 1.10+ apiserver and kubelet.
    - data::Dict{String, String} : Data contains the configuration data. Each key must consist of alphanumeric characters, '-', '_' or '.'. Values with non-UTF-8 byte sequences must use the BinaryData field. The keys stored in Data must not overlap with the keys in the BinaryData field, this is enforced during validation process.
    - kind::String : Kind is a string value representing the REST resource this object represents. Servers may infer this from the endpoint the client submits requests to. Cannot be updated. In CamelCase. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#types-kinds
    - metadata::IoK8sApimachineryPkgApisMetaV1ObjectMeta : Standard object's metadata. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#metadata
"""
mutable struct IoK8sApiCoreV1ConfigMap <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    binaryData::Any # spec type: Union{ Nothing, Dict{String, Vector{UInt8}} } # spec name: binaryData
    data::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: data
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata

    function IoK8sApiCoreV1ConfigMap(;apiVersion=nothing, binaryData=nothing, data=nothing, kind=nothing, metadata=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ConfigMap, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiCoreV1ConfigMap, Symbol("binaryData"), binaryData)
        setfield!(o, Symbol("binaryData"), binaryData)
        validate_property(IoK8sApiCoreV1ConfigMap, Symbol("data"), data)
        setfield!(o, Symbol("data"), data)
        validate_property(IoK8sApiCoreV1ConfigMap, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiCoreV1ConfigMap, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        o
    end
end # type IoK8sApiCoreV1ConfigMap

const _property_map_IoK8sApiCoreV1ConfigMap = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("binaryData")=>Symbol("binaryData"), Symbol("data")=>Symbol("data"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"))
const _property_types_IoK8sApiCoreV1ConfigMap = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("binaryData")=>"Dict{String, Vector{UInt8}}", Symbol("data")=>"Dict{String, String}", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta")
Base.propertynames(::Type{ IoK8sApiCoreV1ConfigMap }) = collect(keys(_property_map_IoK8sApiCoreV1ConfigMap))
Swagger.property_type(::Type{ IoK8sApiCoreV1ConfigMap }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1ConfigMap[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ConfigMap }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ConfigMap[property_name]

function check_required(o::IoK8sApiCoreV1ConfigMap)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ConfigMap }, name::Symbol, val)
end


text:
# reversedimindices (for conversion from row-major to column-major and vice versa)
struct ReverseDimIndices{C<:CartesianIndices}
    revcartinds::C
    length::Int
    ReverseDimIndices(revcartinds::C) where {C<:CartesianIndices} = new{C}(revcartinds, length(revcartinds))
end
ReverseDimIndices(A::AbstractArray) = ReverseDimIndices(CartesianIndices(reverse(size(A))))

function Base.iterate(inds::ReverseDimIndices{N}, state::Int=1) where {N}
    if state > inds.length
        return nothing
    end
    ind = CartesianIndex(reverse(inds.revcartinds[state].I))
    return ind, state + 1
end

Base.length(inds::ReverseDimIndices) = inds.length

reversedimindices(A::AbstractArray) = ReverseDimIndices(A)

# reversedims (for StaticArrays)
reversedims(A::StaticVector) = A
reversedims(A::StaticMatrix) = transpose(A)
# TODO: higher dimensions


text:
@defcomp CH4Cycle begin
    e_globalCH4emissions=Parameter(index=[time],unit="Mtonne/year")
    e_0globalCH4emissions=Parameter(unit="Mtonne/year", default=371.139275895821)
    c_CH4concentration=Variable(index=[time],unit="ppbv")
    pic_preindustconcCH4=Parameter(unit="ppbv", default=700.)
    exc_excessconcCH4=Variable(unit="ppbv")
    c0_CH4concbaseyr=Parameter(unit="ppbv", default=1803.0)
    re_remainCH4=Variable(index=[time],unit="ppbv")
    re_remainCH4base=Variable(unit="ppbv")
    nte_natCH4emissions=Variable(index=[time],unit="Mtonne/year")
    air_CH4fractioninatm=Parameter(unit="%", default=100.)
    tea_CH4emissionstoatm=Variable(index=[time],unit="Mtonne/year")

    # Permafrost outputs (calculations not implemented yet)
    permtce0_permafrostemissions0 = Parameter(unit="Mtonne")
    permtce_permafrostemissions = Parameter(index=[time], unit="Mtonne")

    # Individual teay components
    teaynatural_naturalemissions = Variable(index=[time], unit="Mtonne/t")
    teayanthr_anthropogenicemissions = Variable(index=[time], unit="Mtonne/t")
    teayperm_permafrast = Variable(index=[time], unit="Mtonne/t")

    teay_CH4emissionstoatm=Variable(index=[time],unit="Mtonne/t")
    y_year=Parameter(index=[time],unit="year")
    y_year_0=Parameter(unit="year")
    res_CH4atmlifetime=Parameter(unit="year", default=10.5)
    den_CH4density=Parameter(unit="Mtonne/ppbv", default=2.78)
    stim_CH4emissionfeedback=Parameter(unit="Mtonne/degreeC", default=0.)
    rtl_g0_baselandtemp=Parameter(unit="degreeC", default=1.1683981941310047)
    rtl_g_landtemperature=Parameter(index=[time],unit="degreeC")

    function run_timestep(p, v, d, t)

        if is_first(t)
            #eq.3 from Hope (2006) - natural emissions (carbon cycle) feedback, using global temperatures calculated in ClimateTemperature component
            nte_0=p.stim_CH4emissionfeedback*p.rtl_g0_baselandtemp
            v.nte_natCH4emissions[t]=p.stim_CH4emissionfeedback*p.rtl_g0_baselandtemp
            #eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
            v.tea_CH4emissionstoatm[t]=(p.e_globalCH4emissions[t]+v.nte_natCH4emissions[t])*p.air_CH4fractioninatm/100

            # Individual teay components
            v.teaynatural_naturalemissions[t] = (nte_0+v.nte_natCH4emissions[t])*(p.air_CH4fractioninatm/100)*(p.y_year[t]-p.y_year_0)/2

            tea_0=(p.e_0globalCH4emissions+nte_0)*p.air_CH4fractioninatm/100
            v.teayanthr_anthropogenicemissions[t]=(v.tea_CH4emissionstoatm[t]+tea_0)*(p.y_year[t]-p.y_year_0)/2

            v.teayperm_permafrast[t] = (p.permtce_permafrostemissions[t] - p.permtce0_permafrostemissions0)*(p.air_CH4fractioninatm/100)

            v.teay_CH4emissionstoatm[t] = v.teaynatural_naturalemissions[t] + v.teayanthr_anthropogenicemissions[t] + v.teayperm_permafrast[t]

            #adapted from eq.1 in Hope(2006) - calculate excess concentration in base year
            v.exc_excessconcCH4=p.c0_CH4concbaseyr-p.pic_preindustconcCH4
            #Eq. 2 from Hope (2006) - base-year remaining emissions
            v.re_remainCH4base=v.exc_excessconcCH4*p.den_CH4density
            v.re_remainCH4[t]=v.re_remainCH4base*exp(-(p.y_year[t]-p.y_year_0)/p.res_CH4atmlifetime)+
                v.teay_CH4emissionstoatm[t]*p.res_CH4atmlifetime*(1-exp(-(p.y_year[t]-p.y_year_0)/p.res_CH4atmlifetime))/(p.y_year[t]-p.y_year_0)
        else
            #eq.3 from Hope (2006) - natural emissions (carbon cycle) feedback, using global temperatures calculated in ClimateTemperature component
            #Here assume still using area-weighted average regional temperatures (i.e. land temperatures) for natural emissions feedback
            v.nte_natCH4emissions[t]=p.stim_CH4emissionfeedback*p.rtl_g_landtemperature[t-1] #askChrisHope
            #eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
            v.tea_CH4emissionstoatm[t]=(p.e_globalCH4emissions[t]+v.nte_natCH4emissions[t])*p.air_CH4fractioninatm/100

            # Individual teay components
            v.teaynatural_naturalemissions[t] = (v.nte_natCH4emissions[t-1]+v.nte_natCH4emissions[t])*(p.air_CH4fractioninatm/100)*(p.y_year[t]-p.y_year_0)/2

            #eq.7 from Hope (2006) - average emissions to atm over time period
            v.teayanthr_anthropogenicemissions[t]=(v.tea_CH4emissionstoatm[t]+v.tea_CH4emissionstoatm[t-1])*(p.y_year[t]-p.y_year[t-1])/2

            v.teayperm_permafrast[t] = (p.permtce_permafrostemissions[t] - p.permtce_permafrostemissions[t-1])*(p.air_CH4fractioninatm/100)

            v.teay_CH4emissionstoatm[t] = v.teaynatural_naturalemissions[t] + v.teayanthr_anthropogenicemissions[t] + v.teayperm_permafrast[t]

            #eq.10 from Hope (2006) - remaining emissions in atmosphere
            v.re_remainCH4[t]=v.re_remainCH4[t-1]*exp(-(p.y_year[t]-p.y_year[t-1])/p.res_CH4atmlifetime)+
                v.teay_CH4emissionstoatm[t]*p.res_CH4atmlifetime*(1-exp(-(p.y_year[t]-p.y_year[t-1])/p.res_CH4atmlifetime))/(p.y_year[t]-p.y_year[t-1])
        end

        #eq.11 from Hope(2006) - CH4 concentration
        v.c_CH4concentration[t]=p.pic_preindustconcCH4+v.exc_excessconcCH4*v.re_remainCH4[t]/v.re_remainCH4base
    end
end

function addch4cycle(model::Model, use_permafrost::Bool)
    ch4cycle = add_comp!(model, CH4Cycle)

    if use_permafrost
        ch4cycle[:permtce0_permafrostemissions0] = 934.2010230392067
    else
        ch4cycle[:permtce0_permafrostemissions0] = 0
        ch4cycle[:permtce_permafrostemissions] = zeros(10)
    end

    ch4cycle
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
# SPDX-License-Identifier: MIT

#######################################################################
# Model definition

"""
    struct AffineTransform{D,Q}
        W
        b
    end

Affine transform.
"""
struct AffineTransform{D,Q} <: BMObject
    W::ParameterList
    b::AbstractParameter
end

function (trans::AffineTransform)(x::AbstractVector{<:Number})
    W = hcat([statistics(w) for w in trans.W]...)
    b = statistics(trans.b)
    W' * x + b
end

function (trans::AffineTransform)(X::AbstractVector{<:AbstractVector})
    W = hcat([statistics(w) for w in trans.W]...)
    b = statistics(trans.b)
    [W'*x + b for x in X]
end

#######################################################################
# Hierarchical prior for the bases

function DefaultHNormalDiagParameter(μ::AbstractVector{T}) where T
    EFD.Parameter(μ, identity, identity)
end

struct HNormalDiag{D} <: BMObject
    αs::ParameterList
    param::EFD.AbstractParameter{T} where T
end

function HNormalDiag(αs, μ::AbstractVector{T}) where T
    D = length(μ)
    p = DefaultHNormalDiagParameter(μ)
    HNormalDiag{D}(αs, p)
end

#######################################################################
# EFD.Distribution interface

function _hnormal_split(Tαs)
    ([getindex.(Tαs, 1)...], [getindex.(Tαs, 2)...])
end

function EFD.basemeasure(n::HNormalDiag{D}, x::AbstractVector{T}) where {T,D}
    αs, lnαs = _hnormal_split(statistics.(n.αs))
    sum( T(.5)*(lnαs .- αs .* x.^2 .- log(T(2π))) )
end

function EFD.lognorm(n::HNormalDiag{D},
                     η::AbstractVector{T} = EFD.naturalform(n.param)) where {T,D}
    μ = η
    αs, _ = _hnormal_split(statistics.(n.αs))
    -T(.5)*dot(αs, μ.^2)
end

EFD.stats(::HNormalDiag, x) = x

EFD.splitgrad(n::HNormalDiag, m) = m

function EFD.stdparam(n::HNormalDiag{D},
                      η::AbstractVector{T} = EFD.naturalform(n.param)) where{T,D}
    (μ = η,)
end

function EFD.kldiv(q::EFD.NormalDiag{D}, p::HNormalDiag{D};
                   μ = EFD.gradlognorm(q)) where D
    q_η = EFD.naturalform(q.param)
    αs, lnαs = _hnormal_split(statistics.(p.αs))
    m = EFD.naturalform(p.param)
    T = eltype(m)
    p_η = vcat(αs .* m, T(-.5) .* αs)

    A_p_η = -sum( T(.5) .* (lnαs .- αs .* m.^2 ) )
    A_p_η - EFD.lognorm(q, q_η) - dot(p_η .- q_η, μ)
end

max_stars_repo_name: BUTSpeechFIT/PPCA

max_issues_repo_name: BUTSpeechFIT/PPCA

max_forks_repo_name: BUTSpeechFIT/PPCA


text:
# This example shows how to use your own field implementation for using them to compute on
# elliptic curves. For making things interesting we shall wrap GaloisFields implementeed binary extension field.
using Test
using CryptoGroups
import CryptoGroups: BinaryField, tobits, frombits, specialize
import GaloisFields: ExtensionField, GaloisField

struct GF₂{F <: ExtensionField} <: BinaryField
    x::F
    GF₂{F}(coeffs::NTuple{N, GaloisField(2)}) where {N, F <: ExtensionField} = new(F(coeffs))

    function GF₂{F}(x::BitVector) where F <: ExtensionField
        𝔽₂ = GaloisField(2)
        coeffs = Tuple(𝔽₂(i) for i in x) # May need a reverse
        return GF₂{F}(coeffs)
    end

    GF₂{F}(x::F) where F <: ExtensionField = new{F}(x)
    GF₂(x::F) where F <: ExtensionField = new{F}(x)
end

function specialize(::Type{GF₂}, f::BitVector)
    F, x = GaloisField(2, :x=>f)
    return GF₂{F}
end

specialize(::Type{GF₂}, poly::Vector{Int}) = specialize(GF₂, BitVector(i in poly for i in 0:maximum(poly)))

tobits(a::GF₂) = reverse(BitVector(i.n for i in a.x.coeffs))
frombits(::Type{F}, a::BitVector) where F <: GF₂ = F(reverse(a))

Base.:+(a::F, b::F) where F <: GF₂ = F(a.x + b.x)
Base.:*(a::F, b::F) where F <: GF₂ = F(a.x * b.x)

Base.inv(a::F) where F <: GF₂ = F(inv(a.x))

Base.zero(::Type{GF₂{F}}) where F <: ExtensionField = GF₂(zero(F))
Base.one(::Type{GF₂{F}}) where F <: ExtensionField = GF₂(one(F))

########################## This one we can test easally as follows ##################

import CryptoGroups: Curve_B_163_PB, Curve_K_163_PB, BinaryCurve, generator, oncurve, order, @def, AffinePoint, <|, ECPoint, validate

let
    B_163v3 = specialize(AffinePoint{BinaryCurve, GF₂}, Curve_B_163_PB)

    g = B_163v3 <| generator(Curve_B_163_PB)
    q = order(Curve_B_163_PB)

    @test oncurve(g)
    @test oncurve(g*3)
    @test g * q == zero(g)
end

@def K_163v3 ECPoint{AffinePoint{BinaryCurve, GF₂}} Curve_K_163_PB

let
    g = K_163v3 <| generator(Curve_K_163_PB)

    @test oncurve(g)
    @test oncurve(g*3)
    @test validate(g)
end

max_stars_repo_name: PeaceFounder/DiffieHellmanGroups

max_issues_repo_name: PeaceFounder/DiffieHellmanGroups

max_forks_repo_name: PeaceFounder/DiffieHellmanGroups


text:
using MeshCat
using Colors

using EnhancedGJK
using CoordinateTransformations
using StaticArrays

using GeometryTypes
using GeometryBasics

r1 = GeometryTypes.HyperRectangle(GeometryTypes.Vec(0.0, 0.0, 0.0), GeometryTypes.Vec(1.0, 1.0, 1.0))
r2 = GeometryTypes.HyperRectangle(GeometryTypes.Vec(0.0, 0.0, 0.0), GeometryTypes.Vec(1.0, 1.0, 1.0))

r3 = GeometryBasics.HyperRectangle(GeometryBasics.Vec(0.0, 0.0, 0.0), GeometryBasics.Vec(1.0, 1.0, 1.0))
r4 = GeometryBasics.HyperRectangle(GeometryBasics.Vec(0.0, 0.0, 0.0), GeometryBasics.Vec(1.0, 1.0, 1.0))

s1 = GeometryTypes.HyperSphere(GeometryTypes.Point(0.0, 0.0, 0.0, 0.0), 1.0)
s2 = GeometryBasics.HyperSphere(GeometryBasics.Point(0.0, 0.0, 0.0, 0.0), 1.0)

vis = Visualizer()
open(vis)
p2

result1 = EnhancedGJK.gjk(r1, r2)
result2 = GeometryTypes.gjk(s1, s2)


text:
module AppendedTriangulationsTests

using Test
using Gridap.ReferenceFEs
using Gridap.Arrays
using Gridap.Geometry
using Gridap.Fields
using LinearAlgebra: ⋅

domain = (0,1,0,1)
partition = (10,10)
model = CartesianDiscreteModel(domain,partition)

ncells = num_cells(model)
nin = ceil(Int,2*ncells/3)
cell_to_mask = fill(false,ncells)
cell_to_mask[1:nin] .= true

grid = get_grid(model)

trian_in = Triangulation(model,cell_to_mask)
test_triangulation(trian_in)

trian_out = Triangulation(model,.! cell_to_mask)
test_triangulation(trian_out)

trian = lazy_append(trian_out,trian_in)
test_triangulation(trian)

@test isa(get_cell_coordinates(trian),AppendedArray)

@test isa(get_cell_ref_coordinates(trian),AppendedArray)

@test isa(get_cell_map(trian),AppendedArray)

@test !isa(get_cell_reffe(trian),AppendedArray)

@test isa(get_cell_shapefuns(trian),AppendedArray)

@test isa(get_cell_type(trian),AppendedArray)

glue = get_glue(trian,Val(2))
@test isa(glue.tface_to_mface,AppendedArray)
@test isa(glue.tface_to_mface_map,AppendedArray)

btrian1 = Boundary(model,tags=5)
btrian2 = Boundary(model,tags=2)
btrian = lazy_append(btrian1,btrian2)
test_triangulation(btrian)
@test get_facet_normal(btrian) !== nothing

cell1_to_mat = [ones(3,3) for i in 1:num_cells(btrian1)]
cell2_to_mat = [ones(3,3) for i in 1:num_cells(btrian2)]
cell_to_mat = lazy_append(cell1_to_mat,cell2_to_mat)
tcell_to_mat, ttrian = move_contributions(cell_to_mat,btrian)
@test tcell_to_mat == cell_to_mat
tcell_to_mat, ttrian = move_contributions(collect(cell_to_mat),btrian)
@test tcell_to_mat == cell_to_mat

#order = 1
#quad = CellQuadrature(trian,2*order)
#quad_in = CellQuadrature(trian_in,2*order)
#quad_out = CellQuadrature(trian_out,2*order)
#
#q = get_coordinates(quad)
#w = get_weights(quad)
#@test isa(q,AppendedArray)
#@test isa(w,AppendedArray)
#
## Append triangulations of different cell type
#
#domain = (0,1,0,1)
#partition = (10,10)
#grid1 = CartesianGrid(domain,partition)
#
#domain = (1,2,0,1)
#partition = (10,10)
#grid2 = simplexify(CartesianGrid(domain,partition))
#
#trian = lazy_append(grid1,grid2)
#test_triangulation(trian)

end # module


text:
using Downloads: download
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the linear advection equation

advection_velocity = (0.2, -0.7, 0.5)
equations = LinearScalarAdvectionEquation3D(advection_velocity)

# Solver with polydeg=4 to ensure free stream preservation (FSP) on non-conforming meshes.
# The polydeg of the solver must be at least twice as big as the polydeg of the mesh.
# See https://doi.org/10.1007/s10915-018-00897-9, Section 6.
solver = DGSEM(polydeg=4, surface_flux=flux_lax_friedrichs)

initial_condition = initial_condition_gauss
boundary_condition = BoundaryConditionDirichlet(initial_condition)

boundary_conditions = Dict(
  :all => boundary_condition
)

# Mapping as described in https://arxiv.org/abs/2012.12040 but with less warping.
# The mapping will be interpolated at tree level, and then refined without changing
# the geometry interpolant. The original mapping applied to this unstructured mesh
# causes some Jacobians to be negative, which makes the mesh invalid.
function mapping(xi, eta, zeta)
  # Don't transform input variables between -1 and 1 onto [0,3] to obtain curved boundaries
  # xi = 1.5 * xi_ + 1.5
  # eta = 1.5 * eta_ + 1.5
  # zeta = 1.5 * zeta_ + 1.5

  y = eta + 1/4 * (cos(1.5 * pi * (2 * xi - 3)/3) *
                   cos(0.5 * pi * (2 * eta - 3)/3) *
                   cos(0.5 * pi * (2 * zeta - 3)/3))

  x = xi + 1/4 * (cos(0.5 * pi * (2 * xi - 3)/3) *
                  cos(2 * pi * (2 * y - 3)/3) *
                  cos(0.5 * pi * (2 * zeta - 3)/3))

  z = zeta + 1/4 * (cos(0.5 * pi * (2 * x - 3)/3) *
                    cos(pi * (2 * y - 3)/3) *
                    cos(0.5 * pi * (2 * zeta - 3)/3))

  # Transform the weird deformed cube to be approximately the size of [-5,5]^3 to match IC
  return SVector(5 * x, 5 * y, 5 * z)
end

# Unstructured mesh with 48 cells of the cube domain [-1, 1]^3
mesh_file = joinpath(@__DIR__, "cube_unstructured_2.inp")
isfile(mesh_file) || download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp",
                              mesh_file)

# Mesh polydeg of 2 (half the solver polydeg) to ensure FSP (see above).
mesh = P4estMesh{3}(mesh_file, polydeg=2,
                    mapping=mapping,
                    initial_refinement_level=1)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions=boundary_conditions)

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

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

summary_callback = SummaryCallback()

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

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_restart = SaveRestartCallback(interval=100,
                                   save_final_restart=true)

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

amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first),
                                      base_level=1,
                                      med_level=2, med_threshold=0.1,
                                      max_level=3, max_threshold=0.6)
amr_callback = AMRCallback(semi, amr_controller,
                           interval=5,
                           adapt_initial_condition=true,
                           adapt_initial_condition_only_refine=true)

stepsize_callback = StepsizeCallback(cfl=1.2)

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

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

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


text:
# SPDX-License-Identifier: MIT

struct GSM{T1,T2} <: AbstractModel
    f::Function
    pW::T1
    models::T2
end

function kldiv(gsm::GSM, qθ, Tθ)
    qW, qH = qθ
    μW, μH = Tθ
    kldiv(qW, gsm.pW, μW) + sum(kldiv.(gsm.models, qH, μH))
end

function reparamtrick(x̄, diagx̄x̄ᵀ, vechx̄x̄ᵀ)
    T = eltype(x̄)
    mmᵀ = x̄ * x̄'
    x̄x̄ᵀ = diagm(diagx̄x̄ᵀ) + CompressedSymmetric(size(x̄, 1), 1, vechx̄x̄ᵀ)
    x̄ + cholesky(Symmetric(x̄x̄ᵀ - mmᵀ)).L * randn(T, size(x̄)...)
end

function reparamtrick(x̄::AbstractVector, x̄²::AbstractVector)
    T = eltype(x̄)
    m² = x̄ .^ 2
    x̄ + sqrt.(max.(0, x̄² .- m²)) .* randn(T, size(x̄)...)
end

reparamtrick(μx) = reparamtrick(μx...)

function getstats(gsm::GSM, qθ)
    qW, qH = qθ
    μ(qW), μ.(qH)
end

function unpack(gsm::GSM, qθ, Tθ)
    qW, qH = qθ
    μW, μH = Tθ
    T = eltype(μW)

    #W = reparamtrick(unpack(qW, μW), noisescale=0)
    W = unpack(qW, μW)[1] # Just use the mean, do not learn the variance
    H = hcat(reparamtrick.(unpack.(qH, μH))...)

    Ĥ = vcat(H, ones(T, 1, size(H, 2)))
    gsm.f.(eachcol(W' * Ĥ))
end

loglikelihood(gsm::GSM, Xs, uTθ) = sum(loglikelihood.(gsm.models, Xs, uTθ))

function newposterior(gsm::GSM, qθ::Tuple, ∇μ; lrate=1, clip=true)
    q₀W, q₀H = qθ
    ∇μW, ∇μH = ∇μ

    qW = newposterior(gsm, q₀W, ∇μW; lrate, clip)
    qH = newposterior.(gsm.models, q₀H, ∇μH; lrate, clip)

    qW, qH
end


text:
# test code to check protobuf performance
# compare protobuf serialization performance with julia serialization (benchmark)
#
# last result:
# nloops: 10000
# julia serialization...
# ser byte sz: 38280000
# elapsed time: 10.340241002 seconds (803842820 bytes allocated)
# protobuf serialization...
# ser byte sz: 21670000
# elapsed time: 4.703561939 seconds (433986864 bytes allocated)

module ProtoBufCompareSer
using ProtoBuf
using JSON
using Random
using Serialization

import ProtoBuf.meta

mutable struct TestType
    b::Bool
    i32::Int32
    iu32::UInt32
    i64::Int64
    ui64::UInt64
    f32::Float32
    f64::Float64
    s::AbstractString

    ab::Array{Bool,1}
    ai32::Array{Int32,1}
    ai64::Array{Int64,1}
    af32::Array{Float32,1}
    af64::Array{Float64,1}
    as::Array{AbstractString,1}

    function TestType(fill=false)
        !fill && (return new())
        new(rand(Bool),
            rand(-100:100), rand(1:100),
            rand(-100:100), rand(1:100),
            Float32(rand()*100), Float64(rand()*100),
            randstring(100),
            convert(Array{Bool,1}, rand(Bool,100)),
            round.(Int32, 127*rand(50)),
            round.(Int64, 127*rand(50)),
            rand(Float32, 50),
            rand(Float64, 50),
            [randstring(10) for i in 1:50]
            )
    end
end # type TestType
const __pack_TestType = Symbol[:ab, :ai32, :ai64, :af32, :af64]
meta(t::Type{TestType}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, __pack_TestType)

function julia_ser(t::TestType, n::Int)
    iob = PipeBuffer()
    sz = 0
    for idx in 1:n
        serialize(iob, t)
        sz += iob.size
        deserialize(iob)
    end
    sz
end

function proto_ser(t::TestType, n::Int)
    iob = PipeBuffer()
    sz = 0
    for idx in 1:n
        writeproto(iob, t)
        sz += iob.size
        readproto(iob, TestType())
    end
    sz
end

function json_ser(t::TestType, n::Int)
    iob = PipeBuffer()
    sz = 0
    for idx in 1:n
        JSON.print(iob, t)
        sz += iob.size
        JSON.parse(iob)
    end
    sz
end

println(meta(TestType))

t = TestType(true)
nloops = 10000

println("nloops: $nloops")
GC.gc()
println("julia serialization...")
@time println("ser byte sz: $(julia_ser(t, nloops))")
GC.gc()
println("protobuf serialization...")
@time println("ser byte sz: $(proto_ser(t, nloops))")
GC.gc()
println("JSON serialization...")
@time println("ser byte sz: $(json_ser(t, nloops))")

end # module


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

struct ConceptMapGroup{T <: AbstractConcept}
    owner::T
    concept_maps::Vector{ConceptMap}
end

function ConceptMapGroup(res::Proto.ConceptMapGroup)
    _owner = instantiate(res.owner)
    _concept_maps = ConceptMap(res.concept_maps)
    return ConceptMapGroup(_owner, _concept_maps)
end

function ConceptMapGroup(res::Vector{Proto.ConceptMapGroup})
    _owner = instantiate(res.owner)
    _concept_maps = ConceptMap(res.concept_maps)
    return ConceptMapGroup(_owner, _concept_maps)
end


text:
module tst_cov
    using Test
    using CloudCovErr

    @testset "CovEst" begin
        @test outest_bounds([-1,100],100) == 2
        @test outest_bounds([-1,105],100) == 5

        ref = [
        117  216  315  414  513  612  711  810  909;
        126  225  324  423  522  621  720  819  918;
        135  234  333  432  531  630  729  828  927;
        144  243  342  441  540  639  738  837  936;
        153  252  351  450  549  648  747  846  945;
        162  261  360  459  558  657  756  855  954;
        171  270  369  468  567  666  765  864  963;
        180  279  378  477  576  675  774  873  972;
        189  288  387  486  585  684  783  882  981;
        ]

        arr = reshape(1:121,11,11)
        out = zeros(Int64,9,9)
        tot = zeros(size(arr)[1]);
        boxsmooth!(out, arr, tot, 3, 3)
        @test out == ref

        Np = 5
        widx = 3
        widy = 3
        px0 = 1
        py0 = 1
        sx0 = 8
        sy0 = 8
        tilex = 2
        tiley = 2
        T = Float32
        Δx = (widx-1)÷2
        Δy = (widy-1)÷2
        padx = Np+Δx+px0
        pady = Np+Δy+py0
        stepx = (sx0+2) ÷ tilex
        stepy = (sy0+2) ÷ tiley
        in_subimage = zeros(T,stepx+2*padx,stepy+2*pady)
        in_subimage = convert(Array{Float32,2},reshape(1:size(in_subimage)[1]*size(in_subimage)[2],size(in_subimage)))
        ism = zeros(T,stepx+2*padx,stepy+2*pady)
        bimage = zeros(T,stepx+2*padx-2*Δx,stepy+2*pady-2*Δy)
        bism = zeros(T,stepx+2*padx-2*Δx,stepy+2*pady-2*Δy,2*Np-1, Np);
        cov_avg!(bimage, ism, bism, in_subimage, Np=Np, widx=widx, widy=widy)
        @test bimage[end,end] == 3069
        @test ism[1,end] == 21266.0
        @test bism[:,:,1:4,1] == zeros(stepx+2*padx-2*Δx,stepy+2*pady-2*Δy,4)
        @test bism[end,1,1,2] == 19488.0

        T = Float64
        in_subimage = zeros(T,stepx+2*padx,stepy+2*pady)
        in_subimage = convert(Array{Float32,2},reshape(1:size(in_subimage)[1]*size(in_subimage)[2],size(in_subimage)))
        ism = zeros(T,stepx+2*padx,stepy+2*pady)
        bimage = zeros(T,stepx+2*padx-2*Δx,stepy+2*pady-2*Δy)
        bism = zeros(T,stepx+2*padx-2*Δx,stepy+2*pady-2*Δy,2*Np-1, Np);
        cov_avg!(bimage,ism,bism,in_subimage,Np=Np,widx=widx,widy=widy,ftype=64)
        @test bimage[end,end] == 3069
        @test ism[1,end] == 21266.0
        @test bism[:,:,1:4,1] == zeros(stepx+2*padx-2*Δx,stepy+2*pady-2*Δy,4)
        @test bism[end,1,1,2] == 19488.0
    end
end


text:
"""
k = wavenumber = 2pif/v0
"""
function G0(::T, ::p, ::p, k, m, off::Vararg{Any,3}) where {T<:Acoustic}
    r = sqrt(sum(abs2.(off)))
    G0 = -0.25 * m[:vp] * exp(k * r)
    return G0
end

function G0(::T, ::vz, ::vz, k, m, off::Vararg{Any,3}) where {T<:Elastic}
    r = sqrt(sum(abs2.(off)))
    G0 = -0.25 * m[:vp] * exp(k * r)
    return G0
end
function G0(::T, ::vz, ::vy, k, m, off::Vararg{Any,3}) where {T<:Elastic}
    r = sqrt(sum(abs2.(off)))
    G0 = -0.25 * m[:vp] * exp(k * r)
    return G0
end
function G0(::T, ::vz, ::vx, k, m, off::Vararg{Any,3}) where {T<:Elastic}
    r = sqrt(sum(abs2.(off)))
    G0 = -0.25 * m[:vp] * exp(k * r)
    return G0
end

function G0(::T, ::p, ::p, k, m, off::Vararg{Any,2}) where {T<:Acoustic}
    r = sqrt(sum(abs2.(off)))
    G0 = -0.25 * m[:rho] * im * hankelh2(0, k * r)
    return G0
end

function G0(::T, ::p, ::vz, k, m, off::Vararg{Any,2}) where {T<:Acoustic}
    r = sqrt(sum(abs2.(off)))
    G0 = 0.25 * abs(off[1]) * inv(r) * inv(m[:vp]) * hankelh2(1, k * r)
    return G0
end

function G0(::T, ::p, ::vx, k, m, off::Vararg{Any,2}) where {T<:Acoustic}
    r = sqrt(sum(abs2.(off)))
    G0 = 0.25 * abs(off[2]) * inv(r) * inv(m[:vp]) * hankelh2(1, k * r)
    return G0
end


text:
"""
`block_split(n, max_blk_size)` returns `block_sizes` which divides `n` coordinates into block sizes.

An alternative is to call `block_split` passing
an `n×n` square matrix `Σ` or an `n` vector in place of `n`.
"""
function block_split end

function block_split(n, max_blk_size)
    block_sizes = fill(max_blk_size, n÷max_blk_size)
    m = sum(block_sizes)
    if m < n
        append!(block_sizes, n - m)
    end
    block_sizes
end

function block_split(Σ::AbstractMatrix, max_blk_size)
    n = LinearAlgebra.checksquare(Σ)
    block_split(n, max_blk_size)
end

function block_split(v::AbstractVector, max_blk_size)
    n = length(v)
    block_split(n, max_blk_size)
end


text:
"""
nfetb.jl

NFE-TB Perturbation
"""

immutable NFETB{Tn <: NFEPerturbation, Tt <: TBPerturbation} <: Perturbation
  nfe::Tn # NFE-like perturbation
  tb::Tt # TB-like perturbation
end

ncomp(nfetb::NFETB) = ncomp(nfetb.nfe)
coreradius(nfetb::NFETB) = coreradius(nfetb.nfe)
cutoffradius(nfetb::NFETB) = cutoffradius(nfetb.nfe)
hsdiameter_estimate(nfetb::NFETB) = hsdiameter_estimate(nfetb.nfe)

function pairpotential(nfetb::NFETB)::Array{Function,2}
  u_nfe = pairpotential(nfetb.nfe)
  u_tb = pairpotential(nfetb.tb)

  N = size(u_nfe, 1)
  ret = Array{Function}(N,N)

  for i in 1:N, j in 1:N
    i > j && continue
    u(r) =  u_nfe[i,j](r) + u_tb[i,j](r)
    ret[i,j] = ret[j,i] = u
  end

  return ret
end

function pairpotential_minimizer(nfetb::NFETB)::Array{Float64,2}
  N::Int = ncomp(nfetb)
  rc::Vector{Float64} = nfetb.nfe.pseudo.rc
  u::Array{Function,2} = pairpotential(nfetb)

  ret = Array{Float64,2}(N,N)

  for i in 1:N, j in 1:N
    i > j && continue

    r̄c = (rc[i] + rc[j]) / 2

    opt = Optim.optimize(u[i,j], 3*√r̄c, 7*√r̄c)
    ret[i,j] = ret[j,i] =  Optim.minimizer(opt)
  end

  return ret
end

function pairpotential_derivative(nfetb::NFETB)::Array{Function,2}
  N::Int = ncomp(nfetb)
  u′_nfe::Array{Function,2} = pairpotential_derivative(nfetb.nfe)
  u′_tb::Array{Function,2} = pairpotential_derivative(nfetb.tb)

  [ u′(r) = u′_nfe[i,j](r) + u′_tb[i,j](r) for i in 1:N, j in 1:N ]
end

function entropy(nfetb::NFETB, ref::ReferenceSystem, T::Float64)::Float64
  S_nfe = entropy(nfetb.nfe, ref, T)
  S_tb = entropy(nfetb.tb, ref, T)

  S = S_nfe + S_tb
end

function internal(nfetb::NFETB, ref::ReferenceSystem)::Float64
  U_nfe = internal(nfetb.nfe, ref)
  U_tb = internal(nfetb.tb, ref)

  U = U_nfe + U_tb
end


text:
using OverlapIntegrals, Test

@testset "fact2" begin
    @test fact2(0) == 1
    @test fact2(1) == 1
    @test fact2(2) == 2
    @test fact2(3) == (3 * 1)
    @test fact2(4) == (4 * 2)
    @test fact2(5) == (5 * 3 * 1)
    @test fact2(6) == (6 * 4 * 2)
    @test fact2(7) == (7 * 5 * 3 * 1)
end

@testset "binomial_prefactor" begin
    @test isapprox(binomial_prefactor(1, 1, 1, 0.1, 0.2), 0.3; atol = 1e-7)
    @test isapprox(binomial_prefactor(1, 1, 1, 0.3, 0.4), 0.7; atol = 1e-7)
    @test isapprox(binomial_prefactor(1, 3, 1, 0.1, 0.2), 0.007; atol = 1e-7)
    @test isapprox(binomial_prefactor(2, 3, 1, 0.1, 0.2), 0.09; atol = 1e-7)
end

@testset "overlap1d" begin
    @test isapprox(overlap1d(1, 1, 0.1, 0.2, 1.0), 0.52; atol = 1e-5)
    @test isapprox(overlap1d(3, 1, 0.1, 0.2, 1.0), 0.7952; atol = 1e-5)
end

@testset "tho66" begin
    za = 1.8
    zb = 2.8
    ra = [0.0, 0.0, 0.0]
    rb = [0.5, 0.8, -0.2]

    integral = tho66(za, zb, ra, rb, [0, 0, 0], [0, 0, 0])
    @test isapprox(integral, 0.20373275913014607; atol = 1.0e-16)

    integral = tho66(za, zb, ra, rb, [1, 0, 0], [0, 0, 0])
    @test isapprox(integral, 0.062005622343957505; atol = 1.0e-16)

    integral = tho66(za, zb, ra, rb, [1, 1, 0], [1, 1, 0])
    @test isapprox(integral, -0.00043801221837779696; atol = 1.0e-16)

    integral = tho66(za, zb, ra, rb, [2, 1, 0], [1, 1, 0])
    @test isapprox(integral, -0.0002385994651113168; atol = 1.0e-16)
end

max_stars_repo_name: berquist/overlapintegrals

max_issues_repo_name: berquist/overlapintegrals

max_forks_repo_name: berquist/overlapintegrals


text:
@testset "EarthEnv HabitatHeterogeneity" begin
    using SimpleSDMDataSources: rastername, rasterurl, rasterpath

    @test rastername(EarthEnv{HabitatHeterogeneity}, :cv, 1) == "cv_1km.tif"
    hh_path = joinpath(ENV["ECODATASOURCES_PATH"], "EarthEnv/habitat_heterogeneity")
    @test rasterpath(EarthEnv{HabitatHeterogeneity}) == hh_path
    @test rasterpath(EarthEnv{HabitatHeterogeneity}, :cv, 1) == joinpath(hh_path, "1km/cv_1km.tif")

    @test rasterurl(EarthEnv{HabitatHeterogeneity}, :cv, 1) ==
        "https://data.earthenv.org/habitat_heterogeneity/cv_01_05_1km_uint16.tif"
    download_raster(EarthEnv{HabitatHeterogeneity}; layer=:cv)
    @test isfile(joinpath(hh_path, "25km/cv_25km.tif"))
end


text:
#############################################################################
# sumlargest.jl
# The sum of the k largest/smallest elements of an expression
# All expressions and atoms are subtpyes of AbstractExpr.
# Please read expressions.jl first.
#############################################################################

struct SumLargestAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr}
    size::Tuple{Int, Int}
    k::Int

    function SumLargestAtom(x::AbstractExpr, k::Int)
        if sign(x) == ComplexSign()
            error("Argument should be real instead it is $(sign(x))")
        else
            if k <= 0
                error("sumlargest and sumsmallest only support positive values of k")
            end
            if k > length(x)
                error("k cannot be larger than the number of entries in x")
            end
            children = (x,)
            return new(:sumlargest, hash((children, k)), children, (1,1), k)
        end
    end
end

function sign(x::SumLargestAtom)
    return sign(x.children[1])
end

function monotonicity(x::SumLargestAtom)
    return (Nondecreasing(), )
end

function curvature(x::SumLargestAtom)
    return ConvexVexity()
end

function evaluate(x::SumLargestAtom)
    return sum(sort(vec(evaluate(x.children[1])), rev=true)[1:x.k])
end

function conic_form!(x::SumLargestAtom, unique_conic_forms::UniqueConicForms)
    if !has_conic_form(unique_conic_forms, x)
        c = x.children[1]
        t = Variable(size(c))
        q = Variable()
        # sum k largest given by the solution to
        # minimize sum(t) + k*q
        # subject to c <= t + q, t >= 0
        objective = conic_form!(sum(t) + x.k*q, unique_conic_forms)
        conic_form!(c <= t + q, unique_conic_forms)
        conic_form!(t >= 0, unique_conic_forms)
        cache_conic_form!(unique_conic_forms, x, objective)
    end
    return get_conic_form(unique_conic_forms, x)
end

sumlargest(x::AbstractExpr, k::Int) = k == 0 ? Constant(0) : SumLargestAtom(x, k)
sumsmallest(x::AbstractExpr, k::Int) = k == 0 ? Constant(0) : -SumLargestAtom(-x, k)


text:
type Ray
	origin::Vector
	direction::Vector
end

type Sphere
	center::Vector
	radius::Real
	color::Vector
end

type Intersection
	point::Vector
	distance::Real
	normal::Vector
	object
end


text:
## vec

vec(a::Ones{T}) where T = Ones{T}(length(a))
vec(a::Zeros{T}) where T = Zeros{T}(length(a))
vec(a::Fill{T}) where T = Fill{T}(a.value,length(a))

## Transpose/Adjoint
# cannot do this for vectors since that would destroy scalar dot product

transpose(a::Ones{T,2}) where T = Ones{T}(reverse(a.axes))
adjoint(a::Ones{T,2}) where T = Ones{T}(reverse(a.axes))
transpose(a::Zeros{T,2}) where T = Zeros{T}(reverse(a.axes))
adjoint(a::Zeros{T,2}) where T = Zeros{T}(reverse(a.axes))
transpose(a::Fill{T,2}) where T = Fill{T}(transpose(a.value), reverse(a.axes))
adjoint(a::Fill{T,2}) where T = Fill{T}(adjoint(a.value), reverse(a.axes))

fillsimilar(a::Ones{T}, axes) where T = Ones{T}(axes)
fillsimilar(a::Zeros{T}, axes) where T = Zeros{T}(axes)
fillsimilar(a::AbstractFill, axes) = Fill(getindex_value(a), axes)

permutedims(a::AbstractFill{<:Any,1}) = fillsimilar(a, (1, length(a)))
permutedims(a::AbstractFill{<:Any,2}) = fillsimilar(a, reverse(a.axes))

function permutedims(B::AbstractFill, perm)
    dimsB = size(B)
    ndimsB = length(dimsB)
    (ndimsB == length(perm) && isperm(perm)) || throw(ArgumentError("no valid permutation of dimensions"))
    dimsP = ntuple(i->dimsB[perm[i]], ndimsB)::typeof(dimsB)
    fillsimilar(B, dimsP)
end

## Algebraic identities

function mult_fill(a::AbstractFill, b::AbstractFill{<:Any,2})
    axes(a, 2) ≠ axes(b, 1) &&
        throw(DimensionMismatch("Incompatible matrix multiplication dimensions"))
    return Fill(getindex_value(a)*getindex_value(b), (axes(a, 1), axes(b, 2)))
end

function mult_fill(a::AbstractFill, b::AbstractFill{<:Any,1})
    axes(a, 2) ≠ axes(b, 1) &&
        throw(DimensionMismatch("Incompatible matrix multiplication dimensions"))
    return Fill(getindex_value(a)*getindex_value(b), (axes(a, 1),))
end

function mult_ones(a, b::AbstractMatrix)
    axes(a, 2) ≠ axes(b, 1) &&
        throw(DimensionMismatch("Incompatible matrix multiplication dimensions"))
    return Ones{promote_type(eltype(a), eltype(b))}((axes(a, 1), axes(b, 2)))
end
function mult_ones(a, b::AbstractVector)
    axes(a, 2) ≠ axes(b, 1) &&
        throw(DimensionMismatch("Incompatible matrix multiplication dimensions"))
    return Ones{promote_type(eltype(a), eltype(b))}((axes(a, 1),))
end

function mult_zeros(a, b::AbstractMatrix)
    axes(a, 2) ≠ axes(b, 1) &&
        throw(DimensionMismatch("Incompatible matrix multiplication dimensions"))
    return Zeros{promote_type(eltype(a), eltype(b))}((axes(a, 1), axes(b, 2)))
end
function mult_zeros(a, b::AbstractVector)
    axes(a, 2) ≠ axes(b, 1) &&
        throw(DimensionMismatch("Incompatible matrix multiplication dimensions"))
    return Zeros{promote_type(eltype(a), eltype(b))}((axes(a, 1),))
end

*(a::AbstractFill{<:Any,1}, b::AbstractFill{<:Any,2}) = mult_fill(a,b)
*(a::AbstractFill{<:Any,2}, b::AbstractFill{<:Any,2}) = mult_fill(a,b)
*(a::AbstractFill{<:Any,2}, b::AbstractFill{<:Any,1}) = mult_fill(a,b)

*(a::Ones{<:Any,1}, b::Ones{<:Any,2}) = mult_ones(a, b)
*(a::Ones{<:Any,2}, b::Ones{<:Any,2}) = mult_ones(a, b)
*(a::Ones{<:Any,2}, b::Ones{<:Any,1}) = mult_ones(a, b)

*(a::Zeros{<:Any,1}, b::Zeros{<:Any,2}) = mult_zeros(a, b)
*(a::Zeros{<:Any,2}, b::Zeros{<:Any,2}) = mult_zeros(a, b)
*(a::Zeros{<:Any,2}, b::Zeros{<:Any,1}) = mult_zeros(a, b)

*(a::Zeros{<:Any,1}, b::AbstractFill{<:Any,2}) = mult_zeros(a, b)
*(a::Zeros{<:Any,2}, b::AbstractFill{<:Any,2}) = mult_zeros(a, b)
*(a::Zeros{<:Any,2}, b::AbstractFill{<:Any,1}) = mult_zeros(a, b)
*(a::AbstractFill{<:Any,1}, b::Zeros{<:Any,2}) = mult_zeros(a,b)
*(a::AbstractFill{<:Any,2}, b::Zeros{<:Any,2}) = mult_zeros(a,b)
*(a::AbstractFill{<:Any,2}, b::Zeros{<:Any,1}) = mult_zeros(a,b)

*(a::Zeros{<:Any,1}, b::AbstractMatrix) = mult_zeros(a, b)
*(a::Zeros{<:Any,2}, b::AbstractMatrix) = mult_zeros(a, b)
*(a::AbstractMatrix, b::Zeros{<:Any,1}) = mult_zeros(a, b)
*(a::AbstractMatrix, b::Zeros{<:Any,2}) = mult_zeros(a, b)
*(a::Zeros{<:Any,1}, b::AbstractVector) = mult_zeros(a, b)
*(a::Zeros{<:Any,2}, b::AbstractVector) = mult_zeros(a, b)
*(a::AbstractVector, b::Zeros{<:Any,2}) = mult_zeros(a, b)

*(a::Zeros{<:Any,1}, b::Diagonal) = mult_zeros(a, b)
*(a::Zeros{<:Any,2}, b::Diagonal) = mult_zeros(a, b)
*(a::Diagonal, b::Zeros{<:Any,1}) = mult_zeros(a, b)
*(a::Diagonal, b::Zeros{<:Any,2}) = mult_zeros(a, b)

*(a::Adjoint{T, <:StridedMatrix{T}},   b::Fill{T, 1}) where T = reshape(sum(conj.(parent(a)); dims=1) .* b.value, size(parent(a), 2))
*(a::Transpose{T, <:StridedMatrix{T}}, b::Fill{T, 1}) where T = reshape(sum(parent(a); dims=1) .* b.value, size(parent(a), 2))
*(a::StridedMatrix{T}, b::Fill{T, 1}) where T         = reshape(sum(a; dims=2) .* b.value, size(a, 1))

function *(a::Adjoint{T, <:StridedMatrix{T}}, b::Fill{T, 2}) where T
    fB = similar(parent(a), size(b, 1), size(b, 2))
    fill!(fB, b.value)
    return a*fB
end

function *(a::Transpose{T, <:StridedMatrix{T}}, b::Fill{T, 2}) where T
    fB = similar(parent(a), size(b, 1), size(b, 2))
    fill!(fB, b.value)
    return a*fB
end

function *(a::StridedMatrix{T}, b::Fill{T, 2}) where T
    fB = similar(a, size(b, 1), size(b, 2))
    fill!(fB, b.value)
    return a*fB
end
function *(a::Adjoint{T, <:AbstractVector{T}}, b::Zeros{S, 1}) where {T, S}
    la, lb = length(a), length(b)
    if la ≠ lb
        throw(DimensionMismatch("dot product arguments have lengths $la and $lb"))
    end
    return zero(promote_type(T, S))
end
*(a::Adjoint{T, <:AbstractMatrix{T}} where T, b::Zeros{<:Any, 1}) = mult_zeros(a, b)

function *(a::Transpose{T, <:AbstractVector{T}}, b::Zeros{T, 1}) where T<:Real
    la, lb = length(a), length(b)
    if la ≠ lb
        throw(DimensionMismatch("dot product arguments have lengths $la and $lb"))
    end
    return zero(T)
end
*(a::Transpose{T, <:AbstractMatrix{T}}, b::Zeros{T, 1}) where T<:Real = mult_zeros(a, b)

+(a::Zeros) = a
-(a::Zeros) = a

# Zeros +/- Zeros
function +(a::Zeros{T}, b::Zeros{V}) where {T, V}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    return Zeros{promote_type(T,V)}(size(a)...)
end
-(a::Zeros, b::Zeros) = -(a + b)

# Zeros +/- Fill and Fill +/- Zeros
function +(a::AbstractFill{T}, b::Zeros{V}) where {T, V}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    return convert(AbstractFill{promote_type(T, V)}, a)
end
+(a::Zeros, b::AbstractFill) = b + a
-(a::AbstractFill, b::Zeros) = a + b
-(a::Zeros, b::AbstractFill) = a + (-b)

# Zeros +/- Array and Array +/- Zeros
function +(a::Zeros{T, N}, b::AbstractArray{V, N}) where {T, V, N}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    return AbstractArray{promote_type(T,V),N}(b)
end
function +(a::Array{T, N}, b::Zeros{V, N}) where {T, V, N}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    return AbstractArray{promote_type(T,V),N}(a)
end

function -(a::Zeros{T, N}, b::AbstractArray{V, N}) where {T, V, N}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    return -b + a
end
-(a::Array{T, N}, b::Zeros{V, N}) where {T, V, N} = a + b

+(a::AbstractRange, b::Zeros) = b + a

function +(a::Zeros{T, 1}, b::AbstractRange) where {T}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    Tout = promote_type(T, eltype(b))
    return convert(Tout, first(b)):convert(Tout, step(b)):convert(Tout, last(b))
end
function +(a::Zeros{T, 1}, b::UnitRange) where {T}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    Tout = promote_type(T, eltype(b))
    return convert(Tout, first(b)):convert(Tout, last(b))
end

function -(a::Zeros{T, 1}, b::AbstractRange{V}) where {T, V}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    return -b + a
end
-(a::AbstractRange{T}, b::Zeros{V, 1}) where {T, V} = a + b

####
# norm
####

for op in (:norm1, :norm2, :normInf, :normMinusInf)
    @eval $op(a::Zeros) = norm(getindex_value(a))
end

normp(a::Zeros, p) = norm(getindex_value(a))

norm1(a::AbstractFill) = length(a)*norm(getindex_value(a))
norm2(a::AbstractFill) = sqrt(length(a))*norm(getindex_value(a))
normp(a::AbstractFill, p) = (length(a))^(1/p)*norm(getindex_value(a))
normInf(a::AbstractFill) = norm(getindex_value(a))
normMinusInf(a::AbstractFill) = norm(getindex_value(a))


text:
module Terming

    include("terminal.jl")
    include("event.jl")
    include("parser.jl")
    include("dev_tools.jl")

    export term, in_stream, out_stream, err_stream
    export set_term!

    term = nothing
    in_stream = nothing
    out_stream = nothing
    err_stream = nothing

    function __init__()
        set_term!(init_term())
    end

    function set_term!(t::REPL.Terminals.UnixTerminal)
        global term = t
        global in_stream = t.in_stream
        global out_stream = t.out_stream
        global err_stream = t.err_stream
    end

end


text:
#= polyPlot.jl
Philip Barrett, pobarrett@gmail.com
27may2016, Chicago

Provides plotting functionality for the Polygon class
=#

function polyPlot( pts::Matrix )
  P = [ pts ; pts[1,:] ]
      # Wrap the matrix of points
  if ( size(pts)[1] == 3 ) & (pts[1,1]==pts[2,1]==pts[3,1]) & (pts[1,2]==pts[2,2]==pts[3,2])
    plot( x=P[:,1], y=P[:,2], Geom.point )
  end
  plot( x=P[:,1], y=P[:,2], Geom.path )
end

function polyPlot( poly::Polygon )
  polyPlot( poly.pts )
end

function polyPlot( polys::Array{Polygon,1} )
## TODO: ACCOUNT FOR SINGLETON SETS

  color_vec = ["magenta" "red" "blue" "black" "green" "cyan" "orange"  ]
  # Theme(default_color=)?
  plot( [ layer( x=[ polys[i].pts[:,1] ; polys[i].pts[1,1] ],
                 y=[ polys[i].pts[:,2] ; polys[i].pts[1,2] ],
             Geom.path,
             Theme(default_color=color(parse(Colorant, color_vec[i%7+1]))) )
             for i in 1:length(polys), line_width=8]...)
end


text:
mutable struct HoldActionPolicy{P<:AbstractPolicy} <: AbstractPolicy
    policy::P
    last_action::Vector{Float32}
    action_interval::Int
    counter::Int
end

HoldActionPolicy(; policy::AbstractPolicy, action_interval::Int) =
    HoldActionPolicy(policy, Vector{Float32}(undef, 0), action_interval, 0)

function (p::HoldActionPolicy)(env)
    if p.counter == 0
        p.last_action = p.policy(env)
        p.counter = p.action_interval
    end
    p.counter -= 1
    p.last_action # TODO shoudl this be copied?
end

function RLBase.update!(
    p::HoldActionPolicy,
    traj::CircularArraySARTTrajectory,
    env::AbstractEnv,
    stage::PreActStage,
)
    update!(p.policy, traj, env, stage)
end

log_agent(p::HoldActionPolicy) = log_agent(p.policy)

max_stars_repo_name: albheim/ServiceMeshControl

max_issues_repo_name: albheim/ServiceMeshControl

max_forks_repo_name: albheim/ServiceMeshControl


text:
module Dynamics

using TrajectoryOptimization: Model, UncertainModel, Trajectory, Constraint, Equality, Inequality, Problem
using TrajectoryOptimization
using RigidBodyDynamics
using LinearAlgebra
using DocStringExtensions
using StaticArrays

export
    pendulum,
    doublependulum,
    cartpole,
    cartpole_urdf,
    ballonbeam,
    acrobot_model,
    quadrotor,
    quadrotor_euler,
    kuka,
    doubleintegrator,
    car

include("../dynamics/pendulum.jl")
include("../dynamics/doublependulum.jl")
include("../dynamics/acrobot.jl")
include("../dynamics/ballonbeam.jl")
include("../dynamics/cartpole.jl")
include("../dynamics/quadrotor.jl")
include("../dynamics/quadrotor_euler.jl")
include("../dynamics/kuka.jl")
include("../dynamics/double_integrator.jl")
include("../dynamics/car.jl")

end

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
"""
    struct FermionSpace{N::Int, B <: BasisSet} <: BasisSet

Represents an `N`-particle fermion space, constructed out as a `N`-fold antisymmetrized
tensor product of the single-particle basis set `B`.
"""
struct FermionSpace{N,B <: BasisSet} <: BasisSet
    b :: B

    function FermionSpace(b :: BasisSet, N :: Integer)
        length(b) >= N || error("N=$N too large for basis length(b)=$(length(b)).")
        N = convert(Int, N)
        new{N,typeof(b)}(b)
    end
end
Base.length(fs::FermionSpace{N}) where N = binomial(length(fs.b), N)
nparticles(::FermionSpace{N}) where N = N

struct Fermion1POperator{N, B <: BasisSet, O <: LinearOperator{B}} <: LinearOperator{FermionSpace{N,B}}
    b :: FermionSpace{N,B}
    op :: O
    function Fermion1POperator(fs::FermionSpace{N,B}, op::O) where {N, B <: BasisSet, O <: LinearOperator{B}}
        new{N,B,O}(fs, op)
    end
end
basis(op::Fermion1POperator) = op.b

function apply(i::Integer, op::Fermion1POperator, j::Integer) :: ComplexF64
    N = nparticles(basis(op))
    if i == j
        is = combinationunrank(N, i)
        sum(apply(i, op.op, i) for i = is)
    else
        is, js =  combinationunrank(N, i),  combinationunrank(N, j)
        scd = findcombinationdiff(is, js)
        if isnothing(scd)
            0.0
        else
            idx1, idx2 = scd
            apply(is[idx1], op.op, js[idx2]) * (-1)^(idx1 + idx2)
        end
    end
end


text:
import Term: install_stacktrace

install_stacktrace()

"""
The logic behind these tests is that if something goes
wrong during the error message generation, the
exception returned by the test
will be different from the one you'd expect
"""

@testset "\e[34mERRORS" begin
    @test_throws MethodError 1 - "a"

    @test_throws DomainError √(-1)

    import Term: Panel
    @test_throws AssertionError Panel("mytext", title = "this title is waaaay too long!!!", fit=true)

    @test_throws UndefVarError println(sadfsadfasdsfsd)

    @test_throws BoundsError collect(1:10)[20]

    @test_throws DivideError div(2, 0)

    a() = b()
    b() = a()
    @test_throws StackOverflowError a()

    mydict = Dict(:a => "a", :b => "b")
    @test_throws KeyError mydict["a"]

    @test_throws InexactError Int(2.5)

    function my_func(; my_arg::Int)
        return my_arg + 1
    end
    @test_throws UndefKeywordError my_func()

    m = zeros(20, 20)
    n = zeros(5, 4)
    @test_throws DimensionMismatch m .+ n
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.3.0
#   kernelspec:
#     display_name: Julia 1.3.0
#     language: julia
#     name: julia-1.3
# ---

# + [markdown] slideshow={"slide_type": "slide"}
# # DifferentialEquations.jl
#
# website : https://juliadiffeq.org
#
# *Rackauckas, C. & Nie, Q., (2017). DifferentialEquations.jl – A Performant and Feature-Rich Ecosystem for Solving Differential Equations in Julia. Journal of Open Research Software.*
#
# DifferentialEquations.jl is the core differential equation solver package in Julia. In a nutshell, it provides solvers for:
#
# - ODEs
# - DAEs
# - SODEs, SDAEs
# - Discrete stochastic (Gillespie) equations, mixed with ODEs/SDEs (jump diffusions)
# - DDEs
# - PDEs
#
# It provides wrappers to the common C++ and Fortran libraries, along with native Julia implementations with efficient generic algorithms.
#
# 40 contributors
# -

using DifferentialEquations, Plots

# + [markdown] slideshow={"slide_type": "slide"}
# ## Case Study: Lotka-Volterra in many ways
#
# First, the basic ODE with functions

# + slideshow={"slide_type": "slide"}
function f(du,u,p,t)
  x, y = u
  α, β, γ, δ = p
  du[1] = α*x - β*x*y
  du[2] = -γ*y + δ*x*y
end

p = (1.5,1.0,3.0,1.0)
u0 = [1.0;1.0]
tspan = (0.0,10.0)
prob1 = ODEProblem(f,u0,tspan,p)
sol = solve(prob1)
plot(sol)
# -

plot(sol; vars=(1, 2))

# + [markdown] slideshow={"slide_type": "slide"}
# ## DSLs and Metaprogramming for Syntactic Sugar
#
# Now let's use a DSL defined by a macro:
# -

using ParameterizedFunctions

# +
lorenz63 = @ode_def Lorenz begin
  dx = σ*(y-x)
  dy = ρ*x-y-x*z
  dz = x*y-β*z
end σ β ρ

u0 = [1., 5., 10.]
tspan = (0., 100.)
σ, β, ρ = 10.0, 8.0/3.0, 28.0
p = [σ, β, ρ]
prob = ODEProblem(lorenz63, u0, tspan, p)
sol = solve(prob);
# -

plot(sol, plotdensity=10000, vars=(1, 2, 3))

# + [markdown] slideshow={"slide_type": "slide"}
# ## Optimize the Solver By Choosing from the Largest Set of Methods
#
# http://docs.juliadiffeq.org/latest/solvers/ode_solve.html
#
# https://github.com/JuliaDiffEq/DiffEqBenchmarks.jl
# -

using BenchmarkTools

# + slideshow={"slide_type": "slide"}
@benchmark solve(prob,Vern9())

# + slideshow={"slide_type": "fragment"}
using LSODA
@benchmark solve(prob1,lsoda()) # Default method of SciPy, deSolve, etc.

# + [markdown] slideshow={"slide_type": "slide"}
# ## Generic Algorithms to Propogate Uncertainty

# + slideshow={"slide_type": "fragment"}
using Measurements
u0 = [1.0 ± 0.0 ,1.0 ± 0.0]
p = (1.5 ± 0.0,1.0 ± 0.1,3.0 ± 0.2,1.0 ± 0.1)
tspan = (0.0,10.0)
prob_error = ODEProblem(f,u0,tspan,p)
sol = solve(prob_error,Tsit5(), saveat=0.2)
plot(sol)

# + [markdown] slideshow={"slide_type": "slide"}
# ## Generic Algorithms to Calculate Sensitivities

# + slideshow={"slide_type": "fragment"}
using ForwardDiff: Dual
p1dual = Dual{Float64}(1.5, (1.0, 0.0))
p2dual = Dual{Float64}(1.0, (0.0, 1.0))
pdual = (p1dual, p2dual, 3.0, 1.0)
u0 = [Dual{Float64}(1.0, (0.0, 0.0)),Dual{Float64}(1.0, (0.0, 0.0))]
prob_dual = ODEProblem(f,u0,tspan,pdual)
sol_dual = solve(prob_dual,Tsit5(), saveat=0.2)

timepoints = [i for i in sol_dual.t]
sensitivity_forward_diff = [i[1].partials.values[1] for i in sol_dual.u]
Plots.plot(timepoints,sensitivity_forward_diff,title="dx/da",lw=3,xaxis="Time")

# + [markdown] slideshow={"slide_type": "slide"}
# ## Generic Algorithms for Arbitrary Precision

# + slideshow={"slide_type": "fragment"}
p = big.((1.5,1.0,3.0,1.0))
u0 = big.([1.0,1.0])
tspan = big.((0.0,10.0))
prob1 = ODEProblem(f,u0,tspan,p)
sol = solve(prob1,Vern9(),abstol=1e-25,reltol=1e-25)
sol[10]

# + [markdown] slideshow={"slide_type": "slide"}
# ## Generic Algorithms for Units

# + slideshow={"slide_type": "slide"}
p = (1.5,1.0/u"kg",3.0,1.0/u"kg")./u"s"
u0 = [1.0u"kg",1.0u"kg"]
tspan = (0.0u"s",10.0u"s")
prob_units = ODEProblem(f,u0,tspan,p)
sol = solve(prob_units,Tsit5())
sol.t[10],sol[10]

# + [markdown] slideshow={"slide_type": "slide"}
# # Other Features of DifferentialEquations.jl

# + [markdown] slideshow={"slide_type": "slide"}
# ## The solution is a continuous function by default

# + slideshow={"slide_type": "fragment"}
sol(5.5u"s")

# + [markdown] slideshow={"slide_type": "slide"}
# ## Easily Change to Stochastic Differential Equations

# + slideshow={"slide_type": "fragment"}
function g(du,u,p,t)
  du[1] = 0.2u[1]
  du[2] = 0.2u[2]
end
p = (1.5,1.0,3.0,1.0)
u0 = [1.0;1.0]
tspan = (0.0,10.0)
prob1 = SDEProblem(f,g,u0,tspan,p)
sol = solve(prob1)
plot(sol)

# + [markdown] slideshow={"slide_type": "slide"}
# ## And Delay Differential Equations

# + slideshow={"slide_type": "fragment"}
function f(du,u,h,p,t)
  x, y = u
  α,β,γ,δ = p
  du[1] = α*h(p,t-1)[1] - β*x*y
  du[2] = -γ*y + δ*x*y
end
p = (1.5,1.0,3.0,1.0); u0 = [1.0;1.0]
tspan = (0.0,10.0)
_h(p,t) = ones(2)
prob1 = DDEProblem(f,u0,_h,tspan,p)
sol = solve(prob1)
plot(sol)

max_stars_repo_name: pnavaro/DiffEq

max_issues_repo_name: pnavaro/DiffEq

max_forks_repo_name: pnavaro/DiffEq