text:
using Dates

@testset "Allftplus.jl" begin
    filename = "data\\test1.ALL_FT+"
    filename2 = "data\\test2.ALL_FT+"
    df = DDR2import.Allftplus.read(filename)
    df2 = DDR2import.Allftplus.read(filename2)
    #EXP
    @test df[2,:departureAerodromeIcaoId_0] == "LFST"
    @test df[2,:arrivalAerodromeIcaoId_1] == "EHAM"
    @test df[2,:aircraftId_2] == "AFR27FK"
    @test df[2,:aircraftOperatorIcaoId_3] == "AFR"
    @test df[2,:aircraftTypeIcaoId_4] == "E145"
    @test Dates.year(df[2,:aobt_5]) == 2017
    @test Dates.minute(df[2,:aobt_5]) == 46
    @test df[2,:ifpsId_6] == "AA70843154"
    @test df[2,:aircraftTypeIcaoId_4] == "E145"
    @test Dates.hour(df2[2,:iobt_7]) == 14
    @test df2[2,:originalFlightDataQuality_8] == "FPL"
    @test df2[2,:flightDataQuality_9] == "FPL"
    @test df[2,:source_10] == "UNK"
    @test df[2,:exemptionReasonType_11] == "NEXE"
    @test df[2,:exemptionReasonDistance_12] == "NEXE"
    @test df2[2,:lateFiler_13] == "N"
    @test df2[2,:lateUpdater_14] == "N"
    @test df2[2,:northAtlanticFlight_15] == "N"
    @test Dates.month(df[2,:cobt_16]) == 12
    @test Dates.hour(df[2,:cobt_16]) == 5
    @test Dates.day(df[2,:eobt_17]) == 15
    @test Dates.minute(df[2,:eobt_17]) == 45
    @test Dates.month(df2[2,:lobt_18]) == 1
    @test df[2,:flightState_19] == "NE"
    @test df[2,:previousToActivationFlightState_20] == "NE"
    @test df[2,:suspensionStatus_21] == "NS"
    @test df[2,:tactId_22] == 512256
    @test Dates.year(df[2,:samCtot_23]) == 2017
    @test df2[2,:samSent_24] == "N"
    @test Dates.day(df[2,:sipCtot_25]) == 15
    @test df2[2,:sipSent_26] == "N"
    @test df2[2,:slotForced_27] == "N"
    @test df[2,:mostPenalizingRegulationId_28] == "EHAMA15M"
    @test df[2,:regulationsAffectedByNrOfInstances_29] == 0
    @test df[2,:excludedFromNrOfInstances_30] == 0
    @test df2[6,:lastReceivedAtfmMessageTitle_31] == "REA"
    @test df2[2,:lastReceivedMessageTitle_32] == "FPL"
    @test df2[6,:lastSentAtfmMessageTitle_33] == "SRM"
    @test df[2,:manualExemptionReason_34] == "N"
    @test df2[2,:sensitiveFlight_35] == "N"
    @test df2[2,:readyForImprovement_36] == "Y"
    @test df2[2,:readyToDepart_37] == "N"
    @test df[2,:revisedTaxiTime_38] == 0
    @test df[2,:tis_39] == 0
    @test df[2,:trs_40] == 0
    @test df[2,:toBeSentSlotMessageTitle_41] === missing
    @test df[2,:toBeSentProposalMessageTitle_42] === missing
    @test df2[3,:lastSentSlotMessageTitle_43] === "SRM"
    @test df[2,:lastSentProposalMessageTitle_44] === missing
    @test Dates.year(df[2,:lastSentSlotMessage_45]) == 2017
    @test Dates.day(df[2,:lastSentSlotMessage_45]) == 15
    @test Dates.month(df[2,:lastSentProposalMessage_46]) == 12
    @test df[2,:flightCountOption_47] == "N"
    @test df2[2,:normalFlightTactId_48] == 0
    @test df2[2,:proposalFlightTactId_49] == 0
    @test df2[2,:operatingAircraftOperatorIcaoId_50] == "RSC"
    @test df2[2,:reroutingWhy_51] == "N"
    @test df[2,:reroutedFlightState_52] === missing
    @test df[2,:runwayVisualRange_53] == 0
    @test df2[2,:numberIgnoredErrors_54] == 0
    @test df2[2,:arcAddrSource_55] == "N"
    @test df2[1,:arcAddr_56] == "A9F2C6"
    @test df2[2,:ifpsRegistrationMark_57] == "ECMIF"
    @test df2[2,:flightType_58] == "S"
    @test df2[2,:aircraftEquipment_59] == "DFGRSY:B2B3B4:EB1"
    @test df2[2,:cdmStatus_60] == "N"
    @test Dates.minute(df2[1,:cdmEarlyTtot_61]) == 8
    @test Dates.hour(df2[1,:cdmAoTtot_62]) == 21
    @test Dates.day(df2[1,:cdmAtcTtot_63]) == 2
    @test Dates.month(df2[1,:cdmSequencedTtot_64]) == 1
    @test Dates.minute(df2[1,:cdmTaxiTime_65]) == 21
    @test df2[2,:cdmOffBlockTimeDiscrepancy_66] == "N"
    @test df2[1,:cdmDepartureProcedureId_67] == "DENUT6C"
    @test df2[1,:cdmAircraftTypeId_68] == "B744"
    @test df2[1,:cdmRegistrationMark_69] == "N740CK"
    @test Dates.minute(df2[1,:cdmNoSlotBefore_70]) == 9
    @test df2[2,:cdmDepartureStatus_71] == "K"
    @test df2[2,:ftfmEetFirNrOfInstances_72] == 0
    @test df2[1,:ftfmEetFirList_73] == "EHAAFIR:10"
    @test df2[2,:ftfmEetPtNrOfInstances_74] == 0
    #TODO Setup function to extract EetPtList
    @test df2[1,:ftfmEetPtList_75] == "55N020W:125 57N030W:170 58N040W:211 58N050W:251 CUDDY:281"
    @test df2[2,:ftfmAiracCycleReleaseNumber_76] == 447
    @test df2[2,:ftfmEnvBaselineNumber_77] == 842
    @test df2[2,:ftfmDepartureRunway_78] == "GCXO30"
    @test df2[2,:ftfmArrivalRunway_79] == "GCHI34"
    @test df[2,:ftfmReqFlightlevelSpeedNrOfInstances_80] == 1
    @test df2[1,:ftfmReqFlightlevelSpeedList_81] == "F320:N0487:0 F340:M084:1383 F360:M084:4022"
    fll = DDR2import.Allftplus.reqFlightlevelSpeedList(df2[1,:ftfmReqFlightlevelSpeedList_81])
    @test fll[1].FL == 320
    @test fll[2].Spd == "M084"
    @test fll[3].Value == 4022  # TODO What is Value? Time in seconds? Travel distance in NM?
    @test df2[2,:ftfmConsumedFuel_82] == 322.0
    @test df2[2,:ftfmRouteCharges_83] == 49.0
    @test df2[2,:ftfmAllFtPointNrOfInstances_84] == 16
    pp = DDR2import.Allftplus.AllFtPointProfileList(df2[2,:ftfmAllFtPointProfile_85])
    @test Dates.minute(pp[1].datetime) == 10
    @test pp[5].point == "TESEL"
    @test pp[1].route == "HIE1J"
    @test pp[2].FL == 25
    @test pp[3].pointDistance == 11
    @test pp[4].pointType == "V"
    @test pp[6].geoPointId.ϕ ≈ 28.461666 atol = 0.0001
    @test pp[6].geoPointId.λ ≈ -16.875555555555557 atol = 0.0001
    @test pp[8].ratio == 54
    @test pp[9].isVisible == true
    #@test df[2,:ftfmAllFtPointProfile_85]
    @test df2[2,:ftfmAllFtAirspaceNrOfInstances_86] == 10
    #@test df[2,:ftfmAllFtAirspaceProfile_87]
    ap = DDR2import.Allftplus.AllFtAirspaceProfileList(df2[2,:ftfmAllFtAirspaceProfile_87])
    @test Dates.minute(ap[1].entry_datetime) == 10
    @test ap[2].sector == "GCXOAXO"
    @test Dates.day(ap[3].exit_datetime) == 3
    @test ap[4].fir == "NAS"
    @test ap[5].entry_geoPointId.ϕ ≈ 28.47222222222222 atol = 0.00001
    @test ap[6].exit_geoPointId.λ ≈ -17.3275 atol = 0.00001
    @test ap[7].entry_FL == 80
    @test ap[8].exit_FL == 80
    @test ap[9].entry_pointDistance == 119
    @test ap[10].exit_pointDistance == 185
    @test df2[2,:ftfmAllFtCircleIntersectionsNrOfInstances_88] == 2
    #@test df[2,:ftfmAllFtCircleIntersections_89]
    @test df2[3,:rtfmAiracCycleReleaseNumber_90] == 447
    @test df2[3,:rtfmEnvBaselineNumber_91] == 1498
    @test df2[3,:rtfmDepartureRunway_92] == "EDDL23R"
    @test df2[3,:rtfmArrivalRunway_93] == "LTAC03R"
    @test df2[3,:rtfmReqFlightlevelSpeedNrOfInstances_94] == 3
    @test df2[3,:rtfmReqFlightlevelSpeedList_95] == "F350:N0425:0 F370:N0439:410 F390:N0453:1497"
    @test df2[3,:rtfmConsumedFuel_96] == 7772.0
    @test df2[3,:rtfmRouteCharges_97] == 1153.0
    @test df2[3,:rtfmAllFtPointNrOfInstances_98] == 68
    #@test df[2,:rtfmAllFtPointProfile_99]
    @test df2[3,:rtfmAllFtAirspaceNrOfInstances_100] == 97
    #@test df[2,:rtfmAllFtAirspaceProfile_101]
    @test df2[3,:rtfmAllFtCircleIntersectionsNrOfInstances_102] == 4
    #@test df[2,:rtfmAllFtCircleIntersections_103]
    @test df2[2,:ctfmAiracCycleReleaseNumber_104] == 447
    @test df2[2,:ctfmEnvBaselineNumber_105] == 1270
    @test df2[2,:ctfmDepartureRunway_106] == "GCXO30"
    @test df2[2,:ctfmArrivalRunway_107] == "GCHI34"
    @test df[2,:ctfmReqFlightlevelSpeedNrOfInstances_108] == 1
    @test df2[1,:ctfmReqFlightlevelSpeedList_109] == "F300:N0487:0 F320:N0487:492 F320:M084:1367 F360:M084:4021"
    @test df2[2,:ctfmConsumedFuel_110] == 343.0
    @test df2[2,:ctfmRouteCharges_111] == 49.0
    @test df2[2,:ctfmAllFtPointNrOfInstances_112] == 28
    #@test df[2,:ctfmAllFtPointProfile_113]
    @test df2[2,:ctfmAllFtAirspaceNrOfInstances_114] == 8
    #@test df[2,:ctfmAllFtAirspaceProfile_115]
    @test df2[2,:ctfmAllFtCircleIntersectionsNrOfInstances_116] == 2
    #@test df[2,:ctfmAllFtCircleIntersections_117]
    @test df2[5,:noCPGCPFReason_118] == "X"
    @test Dates.day(df2[2,:scrObt_119]) == 3
    @test df2[2,:scrConsumedFuel_120] == 322.0
    @test df2[2,:scrRouteCharges_121] === 49.0
    @test df2[2,:scrAllFtPointNrOfInstances_122] == 16
    #@test df[2,:scrAllFtPointProfile_123]
    @test df2[2,:scrAllFtAirspaceNrOfInstances_124] == 10
    #@test df[2,:scrAllFtAirspaceProfile_125]
    @test df2[2,:scrAllFtCircleIntersectionsNrOfInstances_126] == 0
    #@test df[2,:scrAllFtCircleIntersections_127]
    @test Dates.minute(df2[3,:srrObt_128]) == 15
    @test df2[3,:srrConsumedFuel_129] == 7813.0
    @test df2[3,:srrRouteCharges_130] == 1113.0
    @test df2[3,:srrAllFtPointNrOfInstances_131] == 74
    #@test df[2,:srrAllFtPointProfile_132]
    @test df2[3,:srrAllFtAirspaceNrOfInstances_133] == 80
    #@test df[2,:srrAllFtAirspaceProfile_134]
    @test df2[3,:srrAllFtCircleIntersectionsNrOfInstances_135] == 4
    #@test df[2,:srrAllFtCircleIntersections_136]
    @test Dates.minute(df2[2,:surObt_137]) == 54
    @test df2[2,:surConsumedFuel_138] == 343.0
    @test df2[2,:surRouteCharges_139] == 49.0
    @test df2[2,:surAllFtPointNrOfInstances_140] == 28
    #@test df[2,:surAllFtPointProfile_141]
    @test df2[2,:surAllFtAirspaceNrOfInstances_142] == 8
    #@test df[2,:surAllFtAirspaceProfile_143]
    @test df2[3,:surAllFtCircleIntersectionsNrOfInstances_144] == 4
    #@test df[2,:surAllFtCircleIntersections_145]
    @test Dates.minute(df2[1,:dctObt_146]) == 50
    @test df2[1,:dctConsumedFuel_147] == 73966.0
    @test df2[1,:dctRouteCharges_148] == 1864.0
    @test df2[1,:dctAllFtPointNrOfInstances_149] == 40
    #@test df[2,:dctAllFtPointProfile_150]
    @test df2[1,:dctAllFtAirspaceNrOfInstances_151] == 79
    #@test df[2,:dctAllFtAirspaceProfile_152]
    @test df2[1,:dctAllFtCircleIntersectionsNrOfInstances_153] == 4
    #@test df[2,:dctAllFtCircleIntersections_154]
    @test df[2,:cpfObt_155] === missing
    @test df[2,:cpfConsumedFuel_156] === missing
    @test df[2,:cpfRouteCharges_157] === missing
    @test df2[2,:cpfAllFtPointNrOfInstances_158] == 0
    #@test df[2,:cpfAllFtPointProfile_159]
    @test df2[2,:cpfAllFtAirspaceNrOfInstances_160] == 0
    #@test df[2,:cpfAllFtAirspaceProfile_161]
    @test df2[2,:cpfAllFtCircleIntersectionsNrOfInstances_162] == 0
    #@test df[2,:cpfAllFtCircleIntersections_163]
    @test df[2,:aircraftidIATA_164] === missing
    @test df[2,:intentionFlight_165] === missing
    @test df[2,:intentionRelatedRouteAssignmentMethod_166] === missing
    @test df[2,:intentionUID_167] === missing
    @test df[2,:intentionEditionDate_168] === missing
    @test df[2,:intentionSource_169] === missing
    @test df[2,:associatedIntentions_170] === missing
    @test df[2,:enrichmentOutput_171] === missing
    @test df2[2,:eventID_172] == "TTE"
    @test Dates.second(df2[1,:eventTime_173]) == 44
    @test df2[2,:flightVersionNr_174] == 44
    @test df[2,:ftfmNrTvProfiles_175] === missing
    # @test df[2,:ftfmTvProfile_176] === missing
    @test df[2,:rtfmNrTvProfiles_177] === missing
    # @test df[2,:rtfmTvProfile_178] === missing
    @test df[2,:ctfmNrTvProfiles_179] === missing
    # @test df[2,:ctfmTvProfile_180] === missing
    #LAST TEST
    @test df[5,:departureAerodromeIcaoId_0] == "WSSS"
end


text:
function sqldf_SQLite(query::String)

    # Normalize Query
    query = replace(query,"\n"=>" ")

    # Extract tables in query
    tables = SQLdf.qtables(query)

    # Remove duplicates
    query = "SELECT * from ("*query*")"

    # Load DataFrame tables in SQLite
    db = SQLite.DB();
    for t in tables
        Main.eval(Main.Meta.parse(t)) |> SQLite.load!(db, t)
    end

    DBInterface.execute(db, query) |> DataFrame

end

ext: jl


text:
	println("	GMTINFO")
	gmt("gmtset -");
	r = gmt("gmtinfo -C", ones(Float32,9,3)*5);
	@assert(r[1].data == [5.0 5 5 5 5 5])
	r = gmtinfo(ones(Float32,9,3)*5, C=true, V=:q);
	@assert(r[1].data == [5.0 5 5 5 5 5])
	#gmtinfo(help=0)

	println("	BLOCK*s")
	d = [0.1 1.5 1; 0.5 1.5 2; 0.9 1.5 3; 0.1 0.5 4; 0.5 0.5 5; 0.9 0.5 6; 1.1 1.5 7; 1.5 1.5 8; 1.9 1.5 9; 1.1 0.5 10; 1.5 0.5 11; 1.9 0.5 12];
	G = blockmedian(region=[0 2 0 2], inc=1, fields="z", reg=true, d);
	if (G !== nothing)	# If run from GMT5 it will return nothing
		G = blockmean(d, region=[0 2 0 2], inc=1, grid=true, reg=true, S=:n);	# Number of points in cell
		G,L = blockmode(region=[0 2 0 2], inc=1, fields="z,l", reg=true, d);
		G,L,H = blockmode(d, region=[0 2 0 2], inc=1, fields="z,l,h", reg=true)
	end
	D = blockmedian(region=[0 2 0 2], inc=1,  reg=true, d);
	D = blockmean(region=[0 2 0 2], inc=1,  reg=true, d);
	D = blockmode(region=[0 2 0 2], inc=1,  reg=true, d);

	println("	CONTOURF")
	G = GMT.peaks();
	gmtwrite("lixo.grd", G)
	C = makecpt(T=(-7,9,2));
	contourf("lixo.grd", Vd=dbg2)
	contourf(G, Vd=dbg2)
	contourf!(G, C=C, Vd=dbg2)
	contourf("lixo.grd", contour=[-2, 0, 2, 5], Vd=dbg2)
	contourf!("lixo.grd", contour=[-2, 0, 2, 5], grd=true, Vd=dbg2)
	contourf(G, C, contour=[-2, 0, 2, 5], Vd=dbg2)
	contourf(G, C, annot=:none, Vd=dbg2)
	d = [0 2 5; 1 4 5; 2 0.5 5; 3 3 9; 4 4.5 5; 4.2 1.2 5; 6 3 1; 8 1 5; 9 4.5 5];
	contourf(d, limits=(-0.5,9.5,0,5), pen=0.25, labels=(line=(:min,:max),), Vd=dbg2)
	C = makecpt(range=(0,10,1));
	contourf(d, C, limits=(-0.5,9.5,0,5), pen=0.25, labels=(line=(:min,:max),), Vd=dbg2)
	contourf(d, C=C, limits=(-0.5,9.5,0,5), pen=0.25, labels=(line=(:min,:max),), Vd=dbg2)

	println("	EARTHTIDE")
	earthtide();
	earthtide(S="");
	earthtide(L=(0,0));

	# FILTER1D
	filter1d([collect((1.0:50)) rand(50)], F="m15");

	println("	FITCIRCLE")
	d = [-3.2488 -1.2735; 7.46259 6.6050; 0.710402 3.0484; 6.6633 4.3121; 12.188 18.570; 8.807 14.397; 17.045 12.865; 19.688 30.128; 31.823 33.685; 39.410 32.460; 48.194 47.114; 62.446 46.528; 59.865 46.453; 68.739 50.164; 64.334 32.984];
	fitcircle(d, L=3);

	println("	GMT2KML & KML2GMT")
	D = gmt("pscoast -R-5/-3/56/58 -Jm1i -M -W0.25p -Di");
	K = gmt2kml(D, F=:l, W=(1,:red));
	gmtwrite("lixo.kml", K)
	kml2gmt("lixo.kml", Z=true);
	kml2gmt(nothing, "lixo.kml", Z=true);	# yes, cheating
	gmtread("lixo.kml", Vd=dbg2);
	#gmtread("lixo.kml");		# Because Travis CI oldness
	rm("lixo.kml")

	println("	GMTCONNECT")
	gmtconnect([0 0; 1 1], [1.1 1.1; 2 2], T=0.5);

	println("	GMTCONVERT")
	gmtconvert([1.1 2; 3 4], o=0)

	println("	GMTGRAVMAG3D")
	if (GMTver > v"6.1.1")
		gmtgravmag3d(M=(shape=:prism, params=(1,1,1,5)), I=1.0, R="-15/15/-15/15", H="10/60/10/-10/40", Vd=dbg2);
		@test_throws ErrorException("Missing one of 'index', 'raw_triang' or 'str' data") gmtgravmag3d(I=1.0);
		@test_throws ErrorException("For grid output MUST specify grid increment ('I' or 'inc')") gmtgravmag3d(Tv=true);
	end

	println("	GMTREGRESS")
	d = [0 5.9 1e3 1; 0.9 5.4 1e3 1.8; 1.8 4.4 5e2 4; 2.6 4.6 8e2 8; 3.3 3.5 2e2 2e1; 4.4 3.7 8e1 2e1; 5.2 2.8 6e1 7e1; 6.1 2.8 2e1 7e1; 6.5 2.4 1.8 1e2; 7.4 1.5 1 5e2];
	regress(d, E=:y, F=:xm, N=2, T="-0.5/8.5/2+n");

	println("	GMTLOGO")
	logo(D="x0/0+w2i")
	logo(julia=8)
	logo(GMTjulia=8, fmt=:png, Vd=dbg2)
	logo(GMTjulia=2, savefig="logo.PNG")
	logo(GMTjulia=2, fmt=:PNG)
	logo!(julia=8, Vd=dbg2)
	logo!("", julia=8, Vd=dbg2)
	@test startswith(logo(pos=(anchor=(0,0),justify=:CM, offset=(1.5,0)), Vd=dbg2), "gmtlogo -Jx1 -Dg0/0+jCM+o1.5/0")
	logo!(julia=8, Vd=dbg2)
	logo!("", julia=8, Vd=dbg2)

	println("	GMTSPATIAL")
	# GMTSPATIAL
	# Test  Cartesian centroid and area
	result = gmt("gmtspatial -Q", [0 0; 1 0; 1 1; 0 1; 0 0]);
	@assert(isapprox(result[1].data, [0.5 0.5 1]))
	# Test Geographic centroid and area
	result = gmt("gmtspatial -Q -fg", [0 0; 1 0; 1 1; 0 1; 0 0]);
	# Intersections
	l1 = gmt("project -C22/49 -E-60/-20 -G10 -Q");
	l2 = gmt("project -C0/-60 -E-60/-30 -G10 -Q");
	#int = gmt("gmtspatial -Ie -Fl", l1, l2);       # Error returned from GMT API: GMT_ONLY_ONE_ALLOWED (59)
	d = [-300 -3500; -200 -800; 400 -780; 500 -3400; -300 -3500];
	gmtspatial(d, C=true, R="0/100/-3100/-3000");

	println("	GMTSELECT")
	gmtselect([2 2], R=(0,3,0,3));		# But is bugged when answer is []
	@test gmtselect([1 2], C=(pts="aa",dist=10), Vd=dbg2) == "gmtselect  -Caa+d10"
	@test gmtselect([1 2], C=(pts=[1 2],dist=10), Vd=dbg2) == "gmtselect  -C+d10"
	@test gmtselect([1 2], C="aa+d0", Vd=dbg2) == "gmtselect  -Caa+d0"
	@test gmtselect([1 2], C=(pts=[1 2],dist=10), L=(line=[1 2;3 4], dist=10), Vd=dbg2) == "gmtselect  -C+d10 -L+d10"

	println("	GMTSET")
	gmtset(MAP_FRAME_WIDTH=0.2)

	println("	GMTSIMPLIFY")
	gmtsimplify([0.0 0; 1.1 1.1; 2 2.2; 3.3 3], T="3k")

	println("	GMTREADWRITE")
	G=gmt("grdmath", "-R0/10/0/10 -I1 5");
	gmtwrite("lixo.grd", G,  scale=10, offset=-10)
	GG = gmtread("lixo.grd", grd=true, varname=:z);
	GG = gmtread("lixo.grd", varname=:z);
	GG = gmtread("lixo.grd", grd=true, layout=:TR);
	GG = gmtread("lixo.grd", grd=true, layout=:TC);
	#GG = gmtread("lixo.grd", grd=true, layer=1);	# This crashes GMT or GDAL in Linux
	@test(sum(G.z[:] - GG.z[:]) == 0)

	GG = mat2grid(reshape(collect(1.0:25), 5, 5));
	gmtwrite("lixo.grd",GG)
	Gr = gmtread("lixo.grd", layout="BCB");
	@test Gr.z[1:5] == [1.0, 2.0, 3.0, 4.0, 5.0]
	Gr = gmtread("lixo.grd", layout="TCB");
	@test Gr.z[1:5] == [5.0, 4.0, 3.0, 2.0, 1.0]
	Gr = gmtread("lixo.grd", layout="TRB");
	@test Gr.z[1:5] == [5.0, 10.0, 15.0, 20.0, 25.0]
	Gr = gmtread("lixo.grd", layout="BRB");
	@test Gr.z[1:5] == [1.0, 6.0, 11.0, 16.0, 21.0]

	gmtwrite("lixo.grd", rand(5,5), id=:cf, layout=:TC)
	gmtwrite("lixo.tif", rand(UInt8,32,32,3), driver=:GTiff)
	I = gmtread("lixo.tif", img=true, layout="TCP");
	I = gmtread("lixo.tif", img=true, layout="ICP");
	I = gmtread("lixo.tif", img=true, band=0);
	I = gmtread("lixo.tif", img=true, band=[0 1 2]);
	show(I);
	imshow(I, show=false)			# Test this one here because we have a GMTimage at hand
	gmtwrite("lixo.tif", mat2img(rand(UInt8,32,32,3)), driver=:GTiff)
	@test GMT.parse_grd_format(Dict(:nan => 0)) == "+n0"
	@test_throws ErrorException("Input data of unknown data type Int64") GMT.gmtwrite(" ", 1)
	@test_throws ErrorException("Number of bands in the 'band' option can only be 1 or 3") GMT.gmtread("", band=[1 2])
	@test_throws ErrorException("Format code MUST have 2 characters and not bla") GMT.parse_grd_format(Dict(:id => "bla"))
	r = rand(UInt8(0):UInt8(10),10,10);	C=makecpt(range=(0,11,1));	I = mat2img(r, cmap=C);
	@test_throws ErrorException("Must select one input data type (grid, image, dataset, cmap or ps)") GG = gmtread("lixo.gr");
	cpt = makecpt(T="-6/8/1");
	gmtwrite("lixo.cpt", cpt)
	cpt = gmtread("lixo.cpt", cpt=true);
	cpt = gmtread("lixo.cpt", cpt=true, Vd=dbg2);
	cpt = gmtread("lixo.cpt");
	gmtwrite("lixo.dat", mat2ds([1 2 10; 3 4 20]))
	gmtwrite("lixo.dat", convert.(UInt8, [1 2 3; 2 3 4]))
	gmtwrite("lixo.dat", [1 2 10; 3 4 20])
	D = gmtread("lixo.dat", i="0,1s10", table=true);
	@test(sum(D[1].data) == 64.0)
	gmtwrite("lixo.dat", D)
	gmt("gmtwrite lixo.cpt", cpt)		# Same but tests other code chunk in gmt_main.jl
	gmt("gmtwrite lixo.dat", D)
	gmt("write lixo.tif=gd:GTiff", mat2img(rand(UInt8,32,32,3)))
	gmt("grdinfo lixo.tif");
	@test_throws ErrorException("First argument cannot be empty. It must contain the file name to write.") gmtwrite("",[1 2]);

	# This gdalwrite screws deeply. Nex GMT errors with non-sensic reason and not reproducible in REPL
	#gdalwrite("lixo1.gmt", fromWKT("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0),(1 1, 1 2, 2 2, 2 1, 1 1))"))
	gmtread("polyg_hole.gmt")

	println("	GMTVECTOR")
	d = [0 0; 0 90; 135 45; -30 -60];
	gmtvector(d, T=:D, S="0/0", f=:g);

	println("	GMTWICH")
	gmtwhich("lixo.dat", C=true)


text:
using SpecialFunctions: besseli, logfactorial, loggamma
import Distributions
using Distributions: ContinuousMultivariateDistribution, AbstractRNG

log_besseli_term(order, m, x) = -logfactorial(m) - loggamma(m + order + 1) + (2*m + order)*log(x/2)

isunitvec(v) = isapprox(sum(x^2 for x in v), 1.)

function logsumexp(arr)
    maxarr = maximum(arr)
    return maxarr + log(sum(exp.(arr .- maxarr)))
end

function _log_besseli(order, x; n_terms=100)
    log_terms = Vector{Float64}(undef, n_terms)
    for m in 0:n_terms-1
        log_terms[m+1] = log_besseli_term(order, m, x)
    end
    return logsumexp(log_terms)
end

const max_lookup_input = 50000
const log_besseli_0_lookup = Vector{Float64}(undef, max_lookup_input) # order 0 (used by plane rotation)
const log_besseli_1_lookup = Vector{Float64}(undef, max_lookup_input) # order 1 (used by 3D rotation)

function log_besseli(order, x; n_terms=Int(floor(x/2)) + 20)
    if x > 500.0 && (order == 0 || order == 1)
        if x > max_lookup_input
            error("log besseli not available for x > $max_loopup_input, got: x = $x")
        end
        floorx = Int(floor(x))
        if order == 0
            if !isassigned(log_besseli_0_lookup, floorx)
                log_besseli_0_lookup[floorx] = _log_besseli(0, floorx; n_terms=n_terms)
            end
            return log_besseli_0_lookup[floorx]
        elseif order == 1
            if !isassigned(log_besseli_1_lookup, floorx)
                log_besseli_1_lookup[floorx] = _log_besseli(1, floorx; n_terms=n_terms)
            end
            return log_besseli_1_lookup[floorx]
        else
            @assert false
        end
    else
        return log(besseli(order, x))
    end
end

function Distributions._vmflck(p, κ)
    T = typeof(κ)
    hp = T(p/2)
    q = hp - 1
    return q * log(κ) - hp * log(2π) - log_besseli(q, κ)
end

import Distributions: VonMisesFisherSampler, _rand!
import LinearAlgebra: qr!, mul!

struct VonMisesFisherSamplerAlternate <: Distributions.Sampleable{Distributions.Multivariate,Distributions.Continuous}
    p::Int          # the dimension
    κ::Float64
    b::Float64
    x0::Float64
    c::Float64
    Q::Matrix{Float64}
end

function Distributions.VonMisesFisherSampler(μ::Vector{Float64}, κ::Float64)
    p = length(μ)
    b = _vmf_bval(p, κ)
    x0 = (1.0 - b) / (1.0 + b)
    c = κ * x0 + (p - 1) * log1p(-abs2(x0))
    Q = _vmf_rotmat(μ)
    VonMisesFisherSamplerAlternate(p, κ, b, x0, c, Q)
end

function Distributions._rand!(rng::AbstractRNG, spl::VonMisesFisherSamplerAlternate,
                x::AbstractVector, t::AbstractVector)
    w = _vmf_genw(rng, spl)
    p = spl.p
    t[1] = w
    s = 0.0
    for i = 2:p
        t[i] = ti = randn(rng)
        s += abs2(ti)
    end

    # normalize t[2:p]
    r = sqrt((1.0 - abs2(w)) / s)
    for i = 2:p
        t[i] *= r
    end

    # rotate
    mul!(x, spl.Q, t)
    return x
end

Distributions._rand!(rng::AbstractRNG, spl::VonMisesFisherSamplerAlternate, x::AbstractVector) =
    _rand!(rng, spl, x, Vector{Float64}(undef, length(x)))

function Distributions._rand!(rng::AbstractRNG, spl::VonMisesFisherSamplerAlternate, x::AbstractMatrix)
    t = Vector{Float64}(undef, size(x, 1))
    for j = 1:size(x, 2)
        _rand!(rng, spl, view(x,:,j), t)
    end
    return x
end

### Core computation

_vmf_bval(p::Int, κ::Real) = (p - 1) / (2.0κ + sqrt(4 * abs2(κ) + abs2(p - 1)))

function _vmf_genw(rng::AbstractRNG, p, b, x0, c, κ)
    # generate the W value -- the key step in simulating vMF
    #
    #   following movMF's document
    #

    r = (p - 1) / 2.0
    betad = Distributions.Beta(r, r)
    z = rand(rng, betad)
    w = (1.0 - (1.0 + b) * z) / (1.0 - (1.0 - b) * z)
    while κ * w + (p - 1) * log(1 - x0 * w) - c < log(rand(rng))
        z = rand(rng, betad)
        w = (1.0 - (1.0 + b) * z) / (1.0 - (1.0 - b) * z)
    end
    return w::Float64
end

_vmf_genw(rng::AbstractRNG, s::VonMisesFisherSamplerAlternate) =
    _vmf_genw(rng, s.p, s.b, s.x0, s.c, s.κ)

function _vmf_rotmat(u::Vector{Float64})
    # construct a rotation matrix Q
    # s.t. Q * [1,0,...,0]^T --> u
    #
    # Strategy: construct a full-rank matrix
    # with first column being u, and then
    # perform QR factorization
    #

    p = length(u)
    A = zeros(p, p)
    copyto!(view(A,:,1), u)

    # let k the be index of entry with max abs
    k = 1
    a = abs(u[1])
    for i = 2:p
        @inbounds ai = abs(u[i])
        if ai > a
            k = i
            a = ai
        end
    end

    # other columns of A will be filled with
    # indicator vectors, except the one
    # that activates the k-th entry
    i = 1
    for j = 2:p
        if i == k
            i += 1
        end
        A[i, j] = 1.0
    end

    # perform QR factorization
    Q = Matrix(qr!(A).Q)
    if dot(view(Q,:,1), u) < 0.0  # the first column was negated
        for i = 1:p
            @inbounds Q[i,1] = -Q[i,1]
        end
    end
    return Q
end


text:
# core.jl

using ImagePhantoms
using ImagePhantoms: rotate2d, coords # helpers
using Unitful: m

using Test: @test, @testset, @test_throws, @inferred

@testset "methods" begin
    ob = Circle(1)
    show(isinteractive() ? stdout : devnull, ob)
    show(isinteractive() ? stdout : devnull, MIME("text/plain"), ob)

    @test ndims(ob) == 2
    @test eltype(ob) == Int64
    @inferred Object(ob)
    @inferred Object(ob; value = 2//3)
    @test Object(ob) == ob
    @test Object(ob ; value = 3).value == 3

    @test phantom([ob]) isa Function
    @test phantom([ob])(0,0) == 1
    @test phantom(zeros(2), zeros(3), [ob]) == ones(2,3)
    @test phantom(zeros(2,3), zeros(2,3), [ob]) == ones(2,3)

    @test radon([ob]) isa Function
    @test radon([ob])(0,0) == 2
    @test radon(zeros(2), zeros(3), [ob]) == 2*ones(2,3)
    @test radon(zeros(2,3), zeros(2,3), [ob]) == 2*ones(2,3)

    @test spectrum([ob]) isa Function
    @test spectrum([ob])(0,0) ≈ π
    @test spectrum(zeros(2), zeros(3), [ob]) ≈ π*ones(2,3)
    @test spectrum(zeros(2,3), zeros(2,3), [ob]) ≈ π*ones(2,3)

end

@testset "helpers" begin
    @test collect(@inferred rotate2d(2, 1, π/2)) ≈ [1,-2]
    @test all((@inferred rotate2d((2, 1), π/2)) .≈ (1,-2))
    @test (@inferred coords(Square(3m), 9m, 6m)) == (3, 2)
end


text:
"""
    SchemeType

A module that only wraps the enum `SchemeTypeEnum` with values `none`, `BFV`, and `CKKS`, which
indicate the type of encryption scheme. `BFV` refers to the Brakerski/Fan-Vercauteren scheme, `CKKS`
refers to the Cheon-Kim-Kim-Song scheme (sometimes also called `HEAAN` in the literature), and
`none` indicates that no encryption should be used.
"""
module SchemeType
@enum SchemeTypeEnum::UInt8 none=0 bfv=1 ckks=2
end

"""
    EncryptionParameters

Stores settings for use by the encryption schemes, most importantly the polynomial modulus, the
coefficient modulus, and the plaintext modulus. An `EncryptionParameters` object is required to
create a `SEALContext` instance.

See also: [`SEALContext`](@ref)
"""
mutable struct EncryptionParameters <: SEALObject
  handle::Ptr{Cvoid}

  function EncryptionParameters(scheme::SchemeType.SchemeTypeEnum=SchemeType.none)
    handleref = Ref{Ptr{Cvoid}}(C_NULL)
    retval = ccall((:EncParams_Create1, libsealc), Clong,
                   (UInt8, Ref{Ptr{Cvoid}}),
                   scheme, handleref)
    @check_return_value retval
    return EncryptionParameters(handleref[])
  end

  function EncryptionParameters(handle::Ptr{Cvoid})
    object = new(handle)
    finalizer(destroy!, object)
    return object
  end
end

function destroy!(object::EncryptionParameters)
  if isallocated(object)
    @check_return_value ccall((:EncParams_Destroy, libsealc), Clong, (Ptr{Cvoid},), object)
    sethandle!(object, C_NULL)
  end

  return nothing
end

function poly_modulus_degree(enc_param::EncryptionParameters)
  degree = Ref{UInt64}(0)
  retval = ccall((:EncParams_GetPolyModulusDegree, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt64}),
                 enc_param, degree)
  @check_return_value retval
  return Int(degree[])
end

function set_poly_modulus_degree!(enc_param::EncryptionParameters, degree)
  retval = ccall((:EncParams_SetPolyModulusDegree, libsealc), Clong,
                 (Ptr{Cvoid}, UInt64),
                 enc_param, degree)
  @check_return_value retval
  return enc_param
end

function set_coeff_modulus!(enc_param::EncryptionParameters, coeff_modulus)
  coeff_modulus_ptrs = Ptr{Cvoid}[gethandle(c) for c in coeff_modulus]
  retval = ccall((:EncParams_SetCoeffModulus, libsealc), Clong,
                 (Ptr{Cvoid}, UInt64, Ref{Ptr{Cvoid}}),
                 enc_param, length(coeff_modulus), coeff_modulus_ptrs)
  @check_return_value retval
  return enc_param
end

function coeff_modulus(enc_param::EncryptionParameters)
  len = Ref{UInt64}(0)

  # First call to obtain length (modulus result pointer is null)
  retval = ccall((:EncParams_GetCoeffModulus, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt64}, Ptr{Cvoid}),
                 enc_param, len, C_NULL)
  @check_return_value retval

  # Second call to obtain modulus
  modulusptrs = Vector{Ptr{Cvoid}}(undef, len[])
  retval = ccall((:EncParams_GetCoeffModulus, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt64}, Ref{Ptr{Cvoid}}),
                 enc_param, len, modulusptrs)
  @check_return_value retval

  modulus = Modulus[Modulus(ptr) for ptr in modulusptrs]
  return modulus
end

function scheme(enc_param::EncryptionParameters)
  scheme = Ref{UInt8}(0)
  retval = ccall((:EncParams_GetScheme, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt8}),
                 enc_param, scheme)
  @check_return_value retval
  return SchemeType.SchemeTypeEnum(scheme[])
end

function plain_modulus(enc_param::EncryptionParameters)
  handleref = Ref{Ptr{Cvoid}}(C_NULL)
  retval = ccall((:EncParams_GetPlainModulus, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{Ptr{Cvoid}}),
                 enc_param, handleref)
  @check_return_value retval
  return Modulus(handleref[], destroy_on_gc=false)
end

function set_plain_modulus!(enc_param::EncryptionParameters, plain_modulus::Modulus)
  retval = ccall((:EncParams_SetPlainModulus1, libsealc), Clong,
                 (Ptr{Cvoid}, Ptr{Cvoid}),
                 enc_param, plain_modulus)
  @check_return_value retval
  return enc_param
end

function set_plain_modulus!(enc_param::EncryptionParameters, plain_modulus::Integer)
  retval = ccall((:EncParams_SetPlainModulus2, libsealc), Clong,
                 (Ptr{Cvoid}, UInt64),
                 enc_param, plain_modulus)
  @check_return_value retval
  return enc_param
end

function parms_id(enc_param::EncryptionParameters)
  parms_id_ = zeros(UInt64, 4)
  retval = ccall((:EncParams_GetParmsId, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt64}),
                 enc_param, parms_id_)
  @check_return_value retval
  return parms_id_
end

function save!(buffer::DenseVector{UInt8}, length::Integer,
               compr_mode::ComprModeType.ComprModeTypeEnum, enc_param::EncryptionParameters)
  out_bytes = Ref{Int64}(0)
  retval = ccall((:EncParams_Save, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt8}, UInt64, UInt8, Ref{Int64}),
                 enc_param, buffer, length, compr_mode, out_bytes)
  @check_return_value retval
  return Int(out_bytes[])
end
function save!(buffer::DenseVector{UInt8}, length::Integer, enc_param::EncryptionParameters)
  return save!(buffer, length, ComprModeType.default, enc_param)
end
function save!(buffer::DenseVector{UInt8}, enc_param::EncryptionParameters)
  return save!(buffer, length(buffer), enc_param)
end

function save_size(compr_mode, enc_param::EncryptionParameters)
  result = Ref{Int64}(0)
  retval = ccall((:EncParams_SaveSize, libsealc), Clong,
                 (Ptr{Cvoid}, UInt8, Ref{Int64}),
                 enc_param, compr_mode, result)
  @check_return_value retval
  return Int(result[])
end
save_size(enc_param::EncryptionParameters) = save_size(ComprModeType.default, enc_param)

function load!(enc_param::EncryptionParameters, buffer::DenseVector{UInt8}, length)
  in_bytes = Ref{Int64}(0)
  retval = ccall((:EncParams_Load, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt8}, UInt64, Ref{Int64}),
                 enc_param, buffer, length, in_bytes)
  @check_return_value retval
  return Int(in_bytes[])
end
load!(enc_param::EncryptionParameters, buffer::DenseVector{UInt8}) = load!(enc_param, buffer,
                                                                           length(buffer))

function Base.:(==)(enc_param1::EncryptionParameters, enc_param2::EncryptionParameters)
  result = Ref{UInt8}(0)
  retval = ccall((:EncParams_Equals, libsealc), Clong,
                 (Ptr{Cvoid}, Ptr{Cvoid}, Ref{UInt8}),
                 enc_param1, enc_param2, result)
  @check_return_value retval
  return Bool(result[])
end


text:
cost(param::AbstractParameter, data::AbstractData) = cost(model(param), param, data)
function cost(model::Model, param, data::AbstractData; idxs=length(model.syms))
    u0 = get_u0(model, param)
    sim = simulate(model, param, u0)
    x = x_values(data)
    y = y_values(data)
    return sum((sim(x; idxs=idxs) .- y) .^ 2)
end


text:
mutable struct SchellingAgent <: AbstractAgent
    id::Int # The identifier number of the agent
    pos::Tuple{Int,Int} # The x, y location of the agent on a 2D grid
    mood::Bool # whether the agent is happy in its position. (true = happy)
    group::Int # The group of the agent,  determines mood as it interacts with neighbors
end

"""
``` julia
schelling(;
    numagents = 320,
    griddims = (20, 20),
    min_to_be_happy = 3
)
```
Same as in [Schelling's segregation model](@ref).
"""
function schelling(; numagents = 320, griddims = (20, 20), min_to_be_happy = 3)
    @assert numagents < prod(griddims)
    space = GridSpace(griddims, periodic = false)
    properties = Dict(:min_to_be_happy => min_to_be_happy)
    model =
        ABM(SchellingAgent, space; properties = properties, scheduler = random_activation)
    ## populate the model with agents, adding equal amount of the two types of agents
    ## at random positions in the model
    for n in 1:numagents
        agent = SchellingAgent(n, (1, 1), false, n < numagents / 2 ? 1 : 2)
        add_agent_single!(agent, model)
    end
    return model, schelling_agent_step!, dummystep
end

function schelling_agent_step!(agent, model)
    agent.mood == true && return # do nothing if already happy
    minhappy = model.min_to_be_happy
    neighbor_positions = nearby_positions(agent, model)
    count_neighbors_same_group = 0
    ## For each neighbor, get group and compare to current agent's group
    ## and increment count_neighbors_same_group as appropriately.
    for neighbor_pos in neighbor_positions
        pos_contents = ids_in_position(neighbor_pos, model)
        ## Skip iteration if the position is empty.
        length(pos_contents) == 0 && continue
        ## Otherwise, get the first agent in the position...
        agent_id = pos_contents[1]
        ## ...and increment count_neighbors_same_group if the neighbor's group is
        ## the same.
        neighbor_agent_group = model[agent_id].group
        if neighbor_agent_group == agent.group
            count_neighbors_same_group += 1
        end
    end
    ## After counting the neighbors, decide whether or not to move the agent.
    ## If count_neighbors_same_group is at least the min_to_be_happy, set the
    ## mood to true. Otherwise, move the agent to a random position.
    if count_neighbors_same_group ≥ minhappy
        agent.mood = true
    else
        move_agent_single!(agent, model)
    end
    return
end


text:
"""
    DynamicMixture

Dummy struct indicating a Dynamic Mixture Model for Time-Series data.
"""
struct DynamicMixture end

"""
    StaticMixture

Dummy struct indicating a Static Mixture Model for Time-Series data.
"""
struct StaticMixture end

"""
    check_dimensions(Y, F_specs, G_specs)

Internal function for checking the problem dimensions and returning them.

Returns observational dimension `n`, number of replicates `nreps`, state
dimensions `p = [p[1], …, p[k]]`, time-window size `T`, number of clusters `k`,
and a vector mapping the `i`-th observation to its indexes in `Y`, in the form
of `UnitRange`'s. Returning a tuple `n, nreps, p, T, k, index_map`.
"""
function check_dimensions(Y::Matrix{RT},
                          F_specs::Vector{Matrix{RT}},
                          G_specs::Vector{Matrix{RT}}) where RT <: Real

    # Number of clusters
    k = size(F_specs, 1)
    if size(G_specs, 1) != k
        throw(DimensionMismatch("F_specs' length does not match G_specs'."))
    end

    # Observational dimensions
    n = size(F_specs[1], 1)
    for j = 2:k
        if size(F_specs[j], 1) != n
            throw(DimensionMismatch("F_specs present inconsistent observational dimensions."))
        end
    end

    # Evolutional dimensions
    p = [size(G_specs[j], 1) for j = 1:k]
    for j = 1:k
        if size(G_specs[j], 2) != p[j]
            throw(DimensionMismatch("Matrices in G_specs should be square."))
        end
        if size(F_specs[j], 2) != p[j]
            throw(DimensionMismatch("F_specs present inconsistent state dimensions."))
        end
    end

    # Number of replicates
    local nreps
    try
        nreps = Int(size(Y, 2) / n)
    catch e
        if isa(e, InexactError)
            throw(DimensionMismatch("Invalid implicit number of replicates"))
        else
            throw(e)
        end
    end

    # Time window length
    T = size(Y, 1)

    # Index mapping
    index_map = [(n * (l - 1) + 1):(n * (l - 1) + n) for l = 1:nreps]

    return n, nreps, p, T, k, index_map
end

"""
    compute_weights(::DynamicMixture, Y, F_specs, G_specs, θ, ϕ, η)

Compute the membership weights of the model, essentially computing the value
from the E-step of the of the Dynamic Mixture Model.
"""
function compute_weights(::DynamicMixture,
                         Y::Matrix{RT},
                         F_specs::Vector{Matrix{RT}},
                         G_specs::Vector{Matrix{RT}},
                         θ::Vector{Matrix{RT}},
                         ϕ::Vector{Vector{RT}},
                         η::Union{Array{RT,3}, Nothing} = nothing) where RT <: Real

    _, nreps, _, T, k, index_map = check_dimensions(Y, F_specs, G_specs)

    if isnothing(η)
        η = ones(RT, T, nreps, k)
    end

    γ = similar(η)

    V = [diagm(1 ./ ϕ[j]) for j = 1:k]

    for t = 1:T
        for i = 1:nreps
            local_index = index_map[i]
            for j = 1:k
                Fⱼ = F_specs[j]
                θⱼ = θ[j]
                Vⱼ = V[j]
                γ[t,i,j] = pdf(MultivariateNormal(Fⱼ * θⱼ[t,:], Vⱼ), Y[t,local_index])
            end
            γ[t,i,:] /= sum(γ[t,i,:])
        end
    end

    return γ
end

"""
    initialize(::DynamicMixture, Y, F_specs, G_specs)

Initializes the model parameters `θ`, `ϕ`, and `η`. Note that to avoid problems
regarding bad relative scales of observational variances, the `ϕ` are
initialized to `ones`.
"""
function initialize(model::Union{DynamicMixture, StaticMixture},
                    Y::Matrix{RT},
                    F_specs::Vector{Matrix{RT}},
                    G_specs::Vector{Matrix{RT}}) where RT <: Real

    n, nreps, p, T, k, index_map = check_dimensions(Y, F_specs, G_specs)

    ϕ = [ones(n) for _ = 1:k]
    θ = Vector{Matrix{RT}}(undef, k)

    # Step 1: Initialize algorithm-specific variables
    centroids = Vector{Int}(undef, 0)
    candidates = collect(1:nreps)
    distances = Matrix{RT}(undef, nreps, 0)

    # Step 2: Pick the first centroid uniformly at random
    push!(centroids, sample(candidates))

    # Step 3: Pick further centroids based on distance
    for _ = 2:k
        # Remove last picked centroid from candidates
        deleteat!(candidates, findfirst(x -> x == centroids[end], candidates))

        # Add to `distances` the distance between every observation and the last
        # picked centroid
        recent_centroid = Y[:,index_map[centroids[end]]]
        local_distances = [sum(abs.(Y[:,index_map[i]] - recent_centroid))  for i = 1:nreps]
        distances = hcat(distances, local_distances)

        # Get each candidate's biggest distance to the centroids
        weights = Weights(maximum(distances, dims=2)[candidates,1])

        # Obtain the new centroid weighted by the its maximum distances
        push!(centroids, sample(candidates, weights))
    end

    # Step 4: Initialize each cluster's parameters based on MAP estimation

    # TODO: Current ordering is assumed to be arbitrary, which is only true if
    # all F_j and G_j are the same. When it isn't adjust all k models for all k
    # centroids and pick the highest likelihood candidate for each model.

    for j = 1:k
        θ[j] = estimate(Y[:,index_map[centroids[j]]], F_specs[j], G_specs[j], 0.7)[1]
    end

    # Step 5: Compute weights

    η = compute_weights(model, Y, F_specs, G_specs, θ, ϕ)

    return θ, ϕ, η
end

"""
    dynamic_mixture(Y, F_specs, G_specs, δ[; numit, mnumit, ϵ, M])

Uses Stochastic Expectation Maximization to obtain MAP estimates for the Dynamic
Mixture of Dynamic Linear Models.

* `Y` is the matrix of observations with `T` rows and `n * nreps` columns;
* `F_specs` is a vector of observational matrices;
* `G_specs` is a vector of evolutional matrices;
* `δ` is the vector of discount factors for each observational unit;
* `numit` is the number of iterations;
* `M` is the number of simulations in the stochastic step of the Stochastic Expectation Maximization.

Returns the estimated `η`, `θ`, `ϕ` and the implied `W`.
"""
function dynamic_mixture(Y::Matrix{RT},
                         F_specs::Vector{Matrix{RT}},
                         G_specs::Vector{Matrix{RT}},
                         δ::Vector{RT};
                         numit::Int = 10,
                         M::Int = 200) where RT <: Real

    _, nobs, _, T, k, _ = check_dimensions(Y, F_specs, G_specs)

    θ, ϕ, η = initialize(DynamicMixture(), Y, F_specs, G_specs)

    c₀ = ones(k)

    mc_observations = Matrix{Int}(undef, T, k)
    mc_estimates = Matrix{RT}(undef, T, k)

    for _ = 1:numit
        # Step 1: Expectation-step

        weights = compute_weights(DynamicMixture(), Y, F_specs, G_specs, θ, ϕ, η)

        # Step 2: Maximize weights for each observation

        for i = 1:nobs
            # Simulate M observations and obtain the mean of all M estimates
            fill!(mc_estimates, 0.)
            for _ = 1:M
                for t = 1:T
                    mc_observations[t,:] = rand(Multinomial(1, η[t,i,:]))
                end
                c = forward_filter(mc_observations, δ[i], c₀)
                mc_estimates += backwards_estimator(c, δ[i]) / M
            end
            η[:,i,:] = mc_estimates
        end

        # Step 3: Maximize the cluster parameters

        for j = 1:k
            θ[j], V = estimate(Y, F_specs[j], G_specs[j], η[:,:,j], 0.7)
            ϕ[j] = 1 ./ diag(V)
        end

    end

    W = [evolutional_covariances(Y, F_specs[j], G_specs[j], Symmetric(diagm(1 ./ ϕ[j])), η[:,:,j], 0.7) for j = 1:k]

    return η, θ, ϕ, W
end


text:
using libm
srand(1)

X = exp(randn(10_000_000))

macro testsum(name,fn,X)
    quote
        function test(X)
            s = 0.0
            for x in X
                s += $fn(x)
            end
            s
        end

        test($X)
        test($X)

        println($name)
        @time test($X)
        println()
    end
end

syslog1p(x::Float64) = Core.Intrinsics.nan_dom_err(ccall((:log1p,:libm),Float64,(Float64,),x), x)
syslog1p(x::Float32) = Core.Intrinsics.nan_dom_err(ccall((:log1pf,:libm),Float32,(Float32,),x), x)

println("Float64")
println("-------------------")
@testsum "Openlibm:" log1p X
@testsum "System:" syslog1p X

#@testsum "Julia fdlibm:" libm.log X
@testsum "Julia tang:" libm.log1p_tang X
# @testsum "muladd:" libm.log_muladd X

println()
println()

X = float32(X)

println("Float32")
println("-------------------")
@testsum "Openlibm:" log1p X
@testsum "System:" syslog1p X

#@testsum "Julia fdlibm:" libm.log X
@testsum "Julia tang:" libm.log1p_tang X


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

using JuLIP, ACE, JuLIP.Potentials, LinearAlgebra
using ACE.Testing: lsq, get_V0
using LinearAlgebra: qr, cond
using Plots

function rand_config(V; rattle = 0.4, nrepeat = 3)
   at = bulk(:W, cubic=true, pbc=false) * nrepeat
   return rattle!(at, rattle)
end

function get_basis(species; N = 3, maxdeg = 10, rcut = 7.0 )
   rcut = 7.0
   basis = ACE.Utils.rpi_basis(; species=species, N = N, r0 = 2.5,
   maxdeg = maxdeg, rcut = rcut,
   rin = rnn(:W) * 0.6,
   constants = false )
   return basis
end

function calc_rmse(train, Vfit)
   sqe = 0.0
   sqy = 0.0
   for (at, E, F) in train
      sqy += (E/length(at))^2 + norm(mat(F)[:])^2
      Efit = energy(Vfit, at)
      Ffit = forces(Vfit, at)
      sqe += ((E - Efit)/length(at))^2 +
             norm( mat(F)[:] - mat(Ffit)[:] )^2
   end
   return sqrt(sqe / sqy)
end

function trainset(Vref, Ntrain; kwargs...)
   train = []
   for n = 1:Ntrain
      at = rand_config(Vref)
      push!(train, (at = at, E = energy(pot, at), F = forces(pot, at)))
   end
   V0 = get_V0(train)
   @show V0
   train = [ (at = at, E = E - energy(V0, at), F = F)
               for (at, E,  F) in train ]
   return train
end

species = :W
pot = EAM("w_eam4.fs")

train = trainset(pot, 100)

basis = get_basis(:W; N = 3, maxdeg = 12, rcut = 6.0)
c, relrmse = lsq(train, basis)

IP = JuLIP.MLIPs.combine(basis, c)

#FORCE SCATTER PLOT ON TEST SET
F_ip = []
F_eam = []

for i in 1:10
   at = rand_config(pot)
   push!(F_ip, vcat(forces(IP, at)...))
   push!(F_eam, vcat(forces(pot, at)...))
end

scatter(vcat(F_ip...), vcat(F_eam...))


text:
using Pkg, LsqFit, Plots, Optim

x = [ 3.08, 4.61, 6.77, 10.15, 12.31, 15.38, 18.77, 22.77, 25.85, 30.15, 36.00, 45.23, 51.69, 61.85, 75.38, 87.08, 110.5]
y = [2.21, 3.59, 6.08, 10.50, 14.09, 19.34, 28.45, 40.33, 50.00, 60.50, 69.89, 80.11, 83.98, 88.95, 93.37, 95.86, 98.07]
a = [0.5, 1, 2, 4]

function oxygen(x,a)
    # a = exp.(a)
    j= 0:(length(a) -1)
    vec(100*sum(j'.*a'.*x.^j',dims=2)./(4*sum(a'.*x.^j',dims=2)))
end
function optim(x,y,a)
    fit = curve_fit(oxygen,x,y,a)
    p = fit.param
    return p
end

function optimize_fn(a)
    sum((y .- oxygen(x,a)).^2)
end
a = ones(5)*(10.)

oxygen(x,a)
optimize_fn(a)

a = rand(5) .* 10
p4_optim = optimize(optimize_fn,a)
round.(p4_optim.minimizer, digits = 2)'
plot!(x,oxygen(x,p4_optim.minimizer),label="prior")
# exp.(p4_optim.minimizer)
p4 = optim(x,y,a)
p4_optimize = optimize(oxygen, )
a = a[1:3]
p3 = optim(x,y,a)

a = a[1:2]
p2 = optim(x,y,a)

plot(x,y,label="data",seriestype=:scatter)
plot!(x,oxygen(x,p4_optim.minimizer),label="prior")
plot!(x,oxygen(x,a),label="prior")
plot!(x,oxygen(x,p4),fillrange=[oxygen(x,conf)],label="optimization 4j")
plot!(x,oxygen(x,p3),label="optimization 3j")
plot!(x,oxygen(x,p2),label="optimization 2j")


text:
"""
    couple_models(aero::Steady, stru::TypicalSection, flap::SimpleFlap)

Create an aerostructural model using a steady aerodynamics model, a
two-degree of freedom typical section model, and a linear steady-state control
surface model.  This model introduces the freestream velocity ``U``, air density
``\\rho``, and control surface deflection ``\\delta`` as additional parameters.
"""
function couple_models(aero::Steady, stru::TypicalSection, flap::SimpleFlap)
    return (aero, stru, flap)
end

"""
    couple_models(aero::QuasiSteady, stru::TypicalSection)

Create an aerostructural model using a quasi-steady aerodynamics model, a
two-degree of freedom typical section model, and a linear steady-state control
surface model.  This model introduces the freestream velocity ``U``, air density
``\\rho``, and control surface deflection ``\\delta`` as additional parameters.
"""
function couple_models(aero::QuasiSteady, stru::TypicalSection, flap::SimpleFlap)
    return (aero, stru, flap)
end

# --- traits --- #

number_of_additional_parameters(::Type{QuasiSteady{0}}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = 3
coupling_inplaceness(::Type{QuasiSteady{0}}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = OutOfPlace()
coupling_rate_jacobian_type(::Type{QuasiSteady{0}}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = Zeros()
coupling_state_jacobian_type(::Type{QuasiSteady{0}}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = Nonlinear()
coupling_parameter_jacobian_type(::Type{QuasiSteady{0}}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = Nonlinear()
coupling_time_gradient_type(::Type{QuasiSteady{0}}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = Zeros()

number_of_additional_parameters(::Type{QuasiSteady{1}}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = 3
coupling_inplaceness(::Type{QuasiSteady{1}}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = OutOfPlace()
coupling_rate_jacobian_type(::Type{QuasiSteady{1}}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = Zeros()
coupling_state_jacobian_type(::Type{QuasiSteady{1}}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = Nonlinear()
coupling_parameter_jacobian_type(::Type{QuasiSteady{1}}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = Nonlinear()
coupling_time_gradient_type(::Type{QuasiSteady{1}}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = Zeros()

number_of_additional_parameters(::Type{QuasiSteady{2}}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = 3
coupling_inplaceness(::Type{QuasiSteady{2}}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = OutOfPlace()
coupling_rate_jacobian_type(::Type{QuasiSteady{2}}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = Linear()
coupling_state_jacobian_type(::Type{QuasiSteady{2}}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = Nonlinear()
coupling_parameter_jacobian_type(::Type{QuasiSteady{2}}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = Nonlinear()
coupling_time_gradient_type(::Type{QuasiSteady{2}}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = Zeros()

# --- methods --- #

function get_coupling_inputs(aero::QuasiSteady{0}, stru::TypicalSection, flap::SimpleFlap,
    dx, x, p, t)
    # extract state variables
    h, θ, hdot, θdot = x
    # extract parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, cnδ, caδ, cmδ, U, ρ, δ = p
    # local freestream velocity components
    u, v = section_steady_velocities(U, θ)
    # calculate aerodynamic loads
    Na, Aa, Ma = quasisteady0_loads(a, b, ρ, a0, α0, cd0, cm0, u, v)
    # loads due to flap deflections
    Nf, Af, Mf = simpleflap_loads(b, u, ρ, cnδ, caδ, cmδ, δ)
    # total loads
    L = Na + Nf
    M = Ma + Mf
    # return inputs
    return SVector(L, M)
end

function get_coupling_inputs(aero::QuasiSteady{1}, stru::TypicalSection, flap::SimpleFlap,
    dx, x, p, t)
    # extract state variables
    h, θ, hdot, θdot = x
    # extract aerodynamic, structural, and aerostructural parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, cnδ, caδ, cmδ, U, ρ, δ = p
    # local freestream velocity components
    u, v, ω = section_velocities(U, θ, hdot, θdot)
    # calculate aerodynamic loads
    Na, Aa, Ma = quasisteady1_loads(a, b, ρ, a0, α0, cd0, cm0, u, v, ω)
    # loads due to flap deflections
    Nf, Af, Mf = simpleflap_loads(b, u, ρ, cnδ, caδ, cmδ, δ)
    # total loads
    L = Na + Nf
    M = Ma + Mf
    # return inputs
    return SVector(L, M)
end

function get_coupling_inputs(aero::QuasiSteady{2}, stru::TypicalSection, flap::SimpleFlap,
    dx, x, p, t)
    # extract rate variables
    dh, dθ, dhdot, dθdot = dx
    # extract state variables
    h, θ, hdot, θdot = x
    # extract aerodynamic, structural, and aerostructural parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, cnδ, caδ, cmδ, U, ρ, δ = p
    # local freestream velocity components
    u, v, ω = section_velocities(U, θ, hdot, θdot)
    udot, vdot, ωdot = section_accelerations(dhdot, dθdot)
    # aerodynamic loads
    Na, Aa, Ma = quasisteady2_loads(a, b, ρ, a0, α0, cd0, cm0, u, v, ω, vdot, ωdot)
    # loads due to flap deflections
    Nf, Af, Mf = simpleflap_loads(b, u, ρ, cnδ, caδ, cmδ, δ)
    # total loads
    L = Na + Nf
    M = Ma + Mf
    # return inputs
    return SVector(L, M)
end

# --- Performance Overloads --- #

function get_coupling_state_jacobian(aero::QuasiSteady{0}, stru::TypicalSection,
    flap::SimpleFlap, dx, x, p, t)
    # extract aerodynamic, structural, and aerostructural parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, cnδ, caδ, cmδ, U, ρ, δ = p
    # return jacobian
    return quasisteady0_section_state_jacobian(a, b, ρ, a0, U)
end

function get_coupling_state_jacobian(aero::QuasiSteady{1}, stru::TypicalSection,
    flap::SimpleFlap, dx, x, p, t)
    # extract aerodynamic, structural, and aerostructural parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, cnδ, caδ, cmδ, U, ρ, δ = p
    # return jacobian
    return quasisteady1_section_state_jacobian(a, b, ρ, a0, U)
end

function get_coupling_state_jacobian(aero::QuasiSteady{2}, stru::TypicalSection,
    flap::SimpleFlap, dx, x, p, t)
    # extract aerodynamic, structural, and aerostructural parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, cnδ, caδ, cmδ, U, ρ, δ = p
    # return jacobian
    return quasisteady2_section_state_jacobian(a, b, ρ, a0, U)
end

# TODO: State rate and parameter jacobians

# --- convenience methods --- #

function set_additional_parameters!(padd, aero::QuasiSteady, stru::TypicalSection,
    flap::SimpleFlap; U, rho, delta)

    padd[1] = U
    padd[2] = rho
    padd[3] = delta

    return padd
end

function separate_additional_parameters(aero::QuasiSteady, stru::TypicalSection,
    flap::SimpleFlap, padd)

    return (U = padd[1], rho = padd[2], delta = padd[3])
end

# --- plotting --- #

@recipe function f(aero::QuasiSteady, stru::TypicalSection, flap::SimpleFlap,
    dx, x, y, p, t)

    framestyle --> :origin
    grid --> false
    xlims --> (-1.0, 1.0)
    ylims --> (-1.5, 1.5)
    label --> @sprintf("t = %6.3f", t)

    # extract state variables
    h, θ, hdot, θdot = x

    # extract parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, cnδ, caδ, cmδ, U, ρ, δ = p

    xplot = [-(0.5 + a*b)*cos(θ),    (0.5 - a*b)*cos(θ)]
    yplot = [ (0.5 + a*b)*sin(θ)-h, -(0.5 - a*b)*sin(θ)-h]

    return xplot, yplot
end


text:
"""
    FunctionalRegret(mf::Tuple{Vararg{<:Function}})

Metric to track the regret of each action using a discrete set of functions.
"""
mutable struct FunctionalRegret{T<:Tuple{Vararg{<:Function}}} <: AbstractMetric
    mf::T
    regret::Vector{Float64}
    function FunctionalRegret(mf::T) where {T<:Tuple{Vararg{<:Function}}}
        return new{T}(mf, Float64[])
    end
end

function (met::FunctionalRegret)(
    X::AbstractMatrix, a::AbstractVector{<:Int}, r::AbstractMatrix
)
    n = size(X, 2)
    Rmean = zeros(n, length(met.mf))
    for a in 1:length(met.mf)
        Rmean[:, a] = mapslices(met.mf[a], X; dims=1)
    end
    Rg = mapslices(row -> maximum(row) .- row, Rmean; dims=2)
    regret = [Rg[i, a[i]] for i in 1:n]
    return met.regret = vcat(met.regret, regret)
end


text:
#
println("Df) Test of the Dielectronic module with ASF from an internally generated initial-, intermediate and final-state multiplets.")
#
setDefaults("print summary: open", "zzz-Dielectronic.sum")
setDefaults("method: continuum, Galerkin")            ## setDefaults("method: continuum, Galerkin") setDefaults("method: continuum, asymptotic Coulomb")
setDefaults("method: normalization, pure sine")       ## setDefaults("method: normalization, pure Coulomb")    setDefaults("method: normalization, pure sine")

if  false
    # K-LL DR into initially He-like tungsten; comparison with Tu et al. (PP, 2016)
    asfSettings   = AsfSettings(AsfSettings(), scField=Basics.HSField())  # not used
    drSettings    = Dielectronic.Settings([E1], [JAC.UseCoulomb, JAC.UseBabushkin], true, PathwaySelection(), -230., 0., 0., CoulombInteraction())
    grid          = Radial.Grid(Radial.Grid(false), rnt = 4.0e-6, h = 5.0e-2, hp = 0.6e-2, rbox = 10.0)
    setDefaults("unit: rate", "1/s");   setDefaults("unit: strength", "cm^2 eV")

    wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=grid, nuclearModel=Nuclear.Model(74.),
                            initialConfigs      = [Configuration("1s^2")],
                            intermediateConfigs = [Configuration("1s 2s^2"),  Configuration("1s 2s 2p"),  Configuration("1s 2p^2")],
                            finalConfigs        = [Configuration("1s^2 2s"),  Configuration("1s^2 2p")],
                            process = Dierec(), processSettings = drSettings )

    wb = perform(wa)

elseif  false
    # K-LL DR into initially Be-like tungsten; comparison with Tu et al. (PP, 2016)
    asfSettings   = AsfSettings(AsfSettings(), scField=Basics.HSField())  # not used
    drSettings    = Dielectronic.Settings([E1], [JAC.UseCoulomb, JAC.UseBabushkin], true, PathwaySelection(), -230., 0., 0., CoulombInteraction())
    grid          = Radial.Grid(Radial.Grid(false), rnt = 4.0e-6, h = 5.0e-2, hp = 0.6e-2, rbox = 10.0)
    setDefaults("unit: rate", "1/s");   setDefaults("unit: strength", "cm^2 eV")

    wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=grid, nuclearModel=Nuclear.Model(74.),
                            initialConfigs      = [Configuration("1s^2 2s^2")],
                            intermediateConfigs = [Configuration("1s 2s^2 2p^2")],
                            finalConfigs        = [Configuration("1s^2 2s^2 2p")],
                            process = Dierec(), processSettings = drSettings )

    wb = perform(wa)

elseif  false
    # K-LL DR into initially O-like tungsten; comparison with Tu et al. (PP, 2016)
    asfSettings   = AsfSettings(AsfSettings(), scField=Basics.HSField())  # not used
    drSettings    = Dielectronic.Settings([E1], [JAC.UseCoulomb, JAC.UseBabushkin], true, PathwaySelection(true, indexTriples=[(1,1,0)]),
                                          -300., 0., 0., CoulombInteraction())
    grid          = Radial.Grid(Radial.Grid(false), rnt = 4.0e-6, h = 5.0e-2, hp = 0.6e-2, rbox = 10.0)
    setDefaults("unit: rate", "1/s");   setDefaults("unit: strength", "cm^2 eV")

    wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=grid, nuclearModel=Nuclear.Model(74.),
                            initialConfigs      = [Configuration("1s^2 2s^2 2p^4")],
                            intermediateConfigs = [Configuration("1s 2s^2 2p^6")],
                            finalConfigs        = [Configuration("1s^2 2s^2 2p^5")],
                            process = Dierec(), processSettings = drSettings )

    wb = perform(wa)

elseif  true
    # K-LL DR into initially He-like carbon; comparison with Xu et al. (PRA, 2016)
    asfSettings   = AsfSettings(AsfSettings(), scField=Basics.HSField())  # not used
    drSettings    = Dielectronic.Settings([E1, M1], [JAC.UseCoulomb, JAC.UseBabushkin], true,
                                          PathwaySelection(true, indexTriples=[(1,1,0), (1,5,0), (1,6,0), (1,10,0), (1,11,0)]),
                                          0., 0., 0., CoulombInteraction())
    grid          = Radial.Grid(Radial.Grid(false), rnt = 4.0e-6, h = 5.0e-2, hp = 0.6e-2, rbox = 12.0)
    setDefaults("unit: rate", "1/s");   setDefaults("unit: strength", "cm^2 eV")

    wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=grid, nuclearModel=Nuclear.Model(6.),
                            initialConfigs      = [Configuration("1s^2")],
                            intermediateConfigs = [Configuration("1s 2s^2"), Configuration("1s 2s 2p"), Configuration("1s 2p^2"), Configuration("1s 2s 3s"),
                                                   Configuration("1s 2s 3p"), Configuration("1s 2s 3d"), Configuration("1s 2p 3s"), Configuration("1s 2p 3p"),
                                                   Configuration("1s 2p 3d")],
                            finalConfigs        = [Configuration("1s^2 2s"), Configuration("1s^2 2p"), Configuration("1s^2 3s"), Configuration("1s^2 3p"),
                                                   Configuration("1s^2 3d"), Configuration("1s 2s 2p"), Configuration("1s 2p^2")],
                            initialAsfSettings  = asfSettings, finalAsfSettings=asfSettings,
                            process = Dierec(), processSettings = drSettings)

    wb = perform(wa)

elseif  false
    # DR if initially Li-like Beryllium; comparison with Mohamed et al. (PRA, 2002) ... however difficult and the identification of resonances and
    # rates does not work at present
    asfSettings   = AsfSettings(AsfSettings(), scField=Basics.HSField())  # not used
    drSettings    = Dielectronic.Settings([E1, M1], [JAC.UseCoulomb, JAC.UseBabushkin], true, PathwaySelection(), 0., 0., 0., CoulombInteraction())
    grid          = Radial.Grid(Radial.Grid(false), rnt = 4.0e-6, h = 5.0e-2, hp = 0.8e-2, rbox = 10.0)
    setDefaults("unit: rate", "1/s")

    wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=grid, nuclearModel=Nuclear.Model(4.),
                            initialConfigs=[Configuration("1s^2 2s")],
                            intermediateConfigs=[Configuration("1s^2 2s^2"),  Configuration("1s^2 2s 2p"),  Configuration("1s^2 2p^2"),
                                                 Configuration("1s^2 2p 3s"), Configuration("1s^2 2p 3p"),  Configuration("1s^2 2p 3d") ],
                            finalConfigs       =[Configuration("1s^2 2s^2"),  Configuration("1s^2 2s 2p"),  Configuration("1s^2 2p^2"),
                                                 Configuration("1s^2 2s 3s"), Configuration("1s^2 2s 3p"),  Configuration("1s^2 2s 3d")],
                            process = Dierec(), processSettings = drSettings )

    wb = perform(wa)

elseif  false
    # DR strength for sodium-like silicon: Comparison with Schmidt et al. (PRA, 2007) ... need to be  adapted
    drSettings    = Dielectronic.Settings([E1, M1], [JAC.UseCoulomb, JAC.UseBabushkin], true,
                                          PathwaySelection(false, indexTriples=[(1,5,0), (1,6,0)]), 0., 0., 0., CoulombInteraction())

grid = Radial.Grid(Radial.Grid(true), rnt = 2.0e-5,h = 5.0e-2, hp = 5.0e-2, NoPoints = 1500)
wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=grid, nuclearModel=Nuclear.Model(14.0),
                        initialConfigs=[Configuration("[Ne] 3s")],
                        intermediateConfigs=[Configuration("[Ne] 3s^2"), Configuration("[Ne] 3s 3p"), Configuration("[Ne] 3s 3d"),
                                             Configuration("[Ne] 3p 4d"), Configuration("[Ne] 3p 4f") ],
                        finalConfigs  =[Configuration("[Ne] 3s^2"),
                                        Configuration("[Ne] 3s 3p"), Configuration("[Ne] 3s 3d"), Configuration("[Ne] 3s 4d"), Configuration("[Ne] 3s 4f") ],
                        process = Dierec(),
                        processSettings=Dielectronic.Settings([E1, M1, M2, E2], [JAC.UseCoulomb, JAC.UseBabushkin], true,
                                                              false, Tuple{Int64,Int64,Int64}[(1,10,1), (1,10,7)], 0., 0., 0., CoulombInteraction())  )

    wb = perform(wa)

elseif  false
    # Need to be adapted !!!
    refConfigs = [Configuration("[Ne] 3s")]
    fromShells = [Shell("3s")]
    toShells   = [Shell("3s"), Shell("3p"), Shell("3d"), Shell("4s"), Shell("4p"), Shell("4d"), Shell("4f"), Shell("5s")]
    intermediateConfigs = Basics.generateConfigurations(refConfigs, fromShells, toShells)
    grid=JAC.Radial.Grid(true)
    #== grid = Radial.Grid(Radial.Grid(true), rnt = 2.0e-5,h = 5.0e-2, hp = 5.0e-2, NoPoints = 2000)
    wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=grid, nuclearModel=Nuclear.Model(74.0, "Fermi"),
                            configs=[Configuration("1s^2 2s^2"), Configuration("1s^2 2s 2p"), Configuration("1s^2 2p^2"), Configuration("1s^2 2s 3s"),
                                    Configuration("1s^2 2s 3p"), Configuration("1s^2 2s 3d"), Configuration("1s^2 2p 3s"), Configuration("1s^2 2p 3p"),
                                    Configuration("1s^2 2p 3d")],
                            ## configs=[Configuration("1s^2 2s"), Configuration("1s^2 2p"),
                            ##          Configuration("1s^2 7s"), Configuration("1s^2 7p"), Configuration("1s^2 7d"), Configuration("1s^2 7f")],
                            asfSettings=AsfSettings(AsfSettings(), scField=Basics.DFSField())  ) ==#

    grid = Radial.Grid(Radial.Grid(true), rnt = 2.0e-5,h = 5.0e-2, hp = 5.0e-2, NoPoints = 1500)
    wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=grid, nuclearModel=Nuclear.Model(74.0),
                            initialConfigs=[Configuration("1s^2 2s"), Configuration("1s^2 2p"),
                                            Configuration("1s^2 7s"), Configuration("1s^2 7p"), Configuration("1s^2 7d"), Configuration("1s^2 7f")],
                            intermediateConfigs=[Configuration("1s^2 2s^2"), Configuration("1s^2 2s 2p"), Configuration("1s^2 2p 7s"), Configuration("1s^2 2p 7p"),
                                                Configuration("1s^2 2p 7d"), Configuration("1s^2 2p 7f")],
                            finalConfigs  =[Configuration("1s^2 2s^2"), Configuration("1s^2 2s 2p"), Configuration("1s^2 2s 7s"), Configuration("1s^2 2s 7p"),
                                            Configuration("1s^2 2s 7d"), Configuration("1s^2 2s 7f")],
                            process = Dierec(),
                            processSettings=Dielectronic.Settings([E1], [JAC.UseCoulomb, JAC.UseBabushkin], true,
                                                                false, Tuple{Int64,Int64,Int64}[(1,10,1), (1,10,7)], 0., 0., 0., CoulombInteraction())  )

    wb = perform(wa)

end
setDefaults("print summary: close", "")


text:
module Projectors
using ..Coords: AbstractCoords

abstract type AbstractProjector end

"""
Coordinate transformation module
"""
struct Projector{Domain, Codomain} where {Domain, Codomain <: AbstractCoords}
    pipeline::Union{Projector, Nothing}
    domain::Type{Domain}
    codomain::Type{Codomain}
end

prj(x) = 2x
inv(x) = x/2

struct Mask <: AbstractProjector

end

end

max_stars_repo_name: schinavro/taps

max_issues_repo_name: schinavro/taps

max_forks_repo_name: schinavro/taps


text:
# function Mpz(pz_pml::SparseMatrixCSC{Float64,Int64}, lambda::Array{Float64,2}, dz::Float64, dt::Float64, ext::Int64, iflag::Int64)
#     lambda = modExpand(lambda, ext, iflag)
#     (m, n) = size(lambda)
#     a1 = 9/8;   a2 = -1/24;
#     c1 = a1/dz; c2 = a2/dz;
#     C3 = zeros(m*n)
#     denum = zeros(m*n)
#     for ix = 1 : n
#         for iz = 1 : m
#             a = 1/dt + pz_pml[iz,ix]/2
#             b = 1/dt - pz_pml[iz,ix]/2
#             denum[(ix-1)*m+iz] = lambda[iz,ix] / a
#             C3[(ix-1)*m+iz]= b / a
#         end
#     end
#     tmp = spzeros(m, m)
#     tmp[1,1] =  c1; tmp[1,2] = c2;
#     tmp[2,1] = -c1; tmp[2,2] = c1; tmp[2,3] = c2;
#     for iz = 3: m-1
#         tmp[iz,iz-1] = -c1; tmp[iz, iz  ] = c1;
#         tmp[iz,iz-2] = -c2; tmp[iz, iz+1] = c2;
#     end
#     tmp[m,m-1] = -c1; tmp[m,m] = c1; tmp[m, m-2] = -c2
#     MpzBpz = spdiagm(C3)
#     MpzBvz = spdiagm(denum) * kron(speye(n), tmp)
#     return MpzBpz, MpzBvz
# end

function Mpz(pz_pml::SparseMatrixCSC{Float64,Int64}, v::Array{Float64,2}, dz::Float64, dt::Float64, ext::Int64, iflag::Int64)
    v = modExpand(v, ext, iflag)
    (m, n) = size(pz_pml)
    a1 = 9/8;   a2 = -1/24;
    c1 = a1/dz; c2 = a2/dz;
    C3 = zeros(m*n)
    denum = zeros(m*n)
    for ix = 1 : n
        for iz = 1 : m
            a = 1/dt + pz_pml[iz,ix]/2
            b = 1/dt - pz_pml[iz,ix]/2
            denum[(ix-1)*m+iz] = (v[iz,ix])^2 / a
            C3[(ix-1)*m+iz]= b / a
        end
    end
    tmp = spzeros(m, m)
    tmp[1,1] =  c1; tmp[1,2] = c2;
    tmp[2,1] = -c1; tmp[2,2] = c1; tmp[2,3] = c2;
    for iz = 3: m-1
        tmp[iz,iz-1] = -c1; tmp[iz, iz  ] = c1;
        tmp[iz,iz-2] = -c2; tmp[iz, iz+1] = c2;
    end
    tmp[m,m-1] = -c1; tmp[m,m] = c1; tmp[m, m-2] = -c2
    MpzBpz = spdiagm(C3)
    MpzBvz = spdiagm(denum) * kron(speye(n), tmp)
    return MpzBpz, MpzBvz
end

function Mpz_back(pz_pml::SparseMatrixCSC{Float64,Int64}, v::Array{Float64,2}, dz::Float64, dt::Float64, ext::Int64, iflag::Int64)
    v = modExpand(v, ext, iflag)
    (m, n) = size(pz_pml)
    a1 = 9/8;   a2 = -1/24;
    c1 = a1/dz; c2 = a2/dz;
    C3 = zeros(m*n)
    denum = zeros(m*n)
    for ix = 1 : n
        for iz = 1 : m
            a = 1/dt - pz_pml[iz,ix]/2
            b = 1/dt + pz_pml[iz,ix]/2
            denum[(ix-1)*m+iz] = (v[iz,ix])^2 / a
            C3[(ix-1)*m+iz]= b / a
        end
    end
    tmp = spzeros(m, m)
    tmp[1,1] =  c1; tmp[1,2] = c2;
    tmp[2,1] = -c1; tmp[2,2] = c1; tmp[2,3] = c2;
    for iz = 3: m-1
        tmp[iz,iz-1] = -c1; tmp[iz, iz  ] = c1;
        tmp[iz,iz-2] = -c2; tmp[iz, iz+1] = c2;
    end
    tmp[m,m-1] = -c1; tmp[m,m] = c1; tmp[m, m-2] = -c2
    MpzBpz = spdiagm(C3)
    MpzBvz = -spdiagm(denum) * kron(speye(n), tmp)
    return MpzBpz, MpzBvz
end

max_stars_repo_name: WneleiGao/AcousticWaveLab

max_issues_repo_name: WneleiGao/AcousticWaveLab

max_forks_repo_name: WneleiGao/AcousticWaveLab


text:
@testset "Double Integrator" begin
    # Test DoubleIntegratorGame
    p = 2
    d = 3
    model = DoubleIntegratorGame(p=p, d=d)
    @test model.n == 2d*p
    @test model.m == d*p
    @test model.p == p
    @test model.ni == 2d*ones(Int,p)
    @test model.mi == d*ones(Int,p)
    @test model.pu[1] == SVector{d,Int}((1,3,5))
    @test model.pu[2] == SVector{d,Int}((2,4,6))

    @test model.px[1] == SVector{2,Int}((1,3))
    @test model.px[2] == SVector{2,Int}((2,4))

    @test model.pz[1] == SVector{2d,Int}((1,3,5,7,9,11))
    @test model.pz[2] == SVector{2d,Int}((2,4,6,8,10,12))

    @test size(model) == (2d*p, d*p, model.pu, p)

    T = Float64
    x = ones(SVector{model.n,T})
    u = 100*ones(SVector{model.m,T})
    @test (@ballocated dynamics($model, $x, $u)) == 0
    # @btime dynamics($model, $x, $u)
end


text:
using Distributions
using Random

"""
    sample(data::DeltaData)

Generate a matrix containing the model inputs for Delta-method sensitivity analysis
with the information in the `data`. The resulting matrix has `N` rows and `D`
columns, where `D` is the number of parameters and `N` is the number of samples.
"""
function sample(data::DeltaData)

    # check for parameters
    if data.params === nothing
        error("Cannot sample to generate model inputs because there are no parameters")
    end

    # constants
    D = length(data.params) # number of uncertain parameters in problem
    N = data.N # number of samples
    d = 1.0 / N

    # preallocate
    result = zeros(N, D)
    temp = zeros(N)

    # generate samples
    for i = 1:D
        for j = 1:N

            low = (j - 1) * d
            high = j * d
            uniform_distrib = Uniform(low, high)

            temp[j] = rand(uniform_distrib)
        end
        shuffle!(temp)
        result[:, i] = temp
    end

    scale_samples!(result, [values(data.params)...]) # scale
    return result

end

max_stars_repo_name: lrennels/SAJulia

max_issues_repo_name: lrennels/SAJulia


text:
using LocallyWeightedRegression
using GeoStatsBase
using Variography
using Plots; gr()
using VisualRegressionTests
using Test, Pkg, Random

ENV["GKSwstype"] = "100"

# environment settings
islinux = Sys.islinux()
istravis = "TRAVIS" ∈ keys(ENV)
datadir = joinpath(@__DIR__,"data")
visualtests = !istravis || (istravis && islinux)
if !istravis
  Pkg.add("Gtk")
  using Gtk
end

@testset "1D regression problem" begin
  Random.seed!(2017)

  N = 100
  x = range(0, stop=1, length=N)
  y = x.^2 .+ [i/1000*randn() for i=1:N]

  sdata   = PointSetData(Dict(:y => y), reshape(x, 1, length(x)))
  sdomain = RegularGrid((0.,), (1.,), dims=(N,))
  problem = EstimationProblem(sdata, sdomain, :y)

  solver = LocalWeightRegress(:y => (variogram=ExponentialVariogram(range=3/10),))

  solution = solve(problem, solver)

  yhat, yvar = solution[:y]

  if visualtests
    @plottest begin
      scatter(x, y, label="data", size=(1000,400))
      plot!(x, yhat, ribbon=yvar, fillalpha=.5, label="LWR")
    end joinpath(datadir,"solution.png") !istravis
  end
end


text:
"""
    function chain(X, arr::Array{L,1}) where {L<:Layer}

Returns the input `Layer` and the output `Layer` from an `Array` of layers and the input of the model as and `Array` `X`

"""
function chain(X, arr::Array{L,1}) where {L<:Layer}
    if ! (arr[1] isa Input)
        X_Input = Input(X)
    end

    for l=1:length(arr)
        if l==1
            if ! isa(arr[l],Input)
                X = arr[l](X_Input)
            else
                X_Input = X = arr[l](X)
            end
        else
            X = arr[l](X)
        end
    end #for

    return X_Input, X
end

export chain

"""
    connect with the previous layer
"""
function (l::FCLayer)(li_1::Layer)
    l.prevLayer = li_1
    if ! in(l,li_1.nextLayers)
        push!(li_1.nextLayers, l)
    end
    l.inputS = li_1.outputS
    l.outputS = (l.channels, li_1.outputS[2])
    return l
end #function (l::FCLayer)(li_1::Layer)

function (l::Activation)(li_1::Layer)
    l.prevLayer = li_1

    l.channels = li_1.channels

    l.inputS = l.outputS = li_1.outputS

    if ! in(l,li_1.nextLayers)
        push!(li_1.nextLayers, l)
    end
    return l
end

function (l::Flatten)(li_1::Layer)
    l.prevLayer = li_1

    l.inputS = li_1.outputS

    l.channels = prod(l.inputS[1:end-1])

    l.outputS = (l.channels, l.inputS[end])

    if !(l in li_1.nextLayers)
        push!(li_1.nextLayers, l)
    end

    return l
end

function (l::BatchNorm)(li_1::Layer)
    l.prevLayer = li_1

    l.channels = li_1.channels

    N = length(li_1.outputS)

    if l.dim >= N
        throw(DimensionMismatch("Normalization Dimension must be less than $N"))
    end

    l.inputS = l.outputS = li_1.outputS

    if ! in(l,li_1.nextLayers)
        push!(li_1.nextLayers, l)
    end
    return l

end #function (l::BatchNorm)

function (l::ConcatLayer)(ls::Array{L,1}) where {L<:Layer}
    for li in ls
        if !in(li,l.prevLayer)
            push!(l.prevLayer, li)
        end
        if !in(l, li.nextLayers)
            push!(li.nextLayers, l)
        end
    end #for

    l.inputS = ls[1].outputS
    l.channels = sum(li.channels for li in ls)
    l.LSlice[l.prevLayer[1]] = 1:l.prevLayer[1].channels
    for i=2:length(l.prevLayer)
        l.LSlice[l.prevLayer[i]] = (l.prevLayer[i-1].channels+1) : (l.prevLayer[i-1].channels+l.prevLayer[i].channels)
    end
    l.outputS = (l.inputS[1:end-2]...,l.channels,l.inputS[end])

    return l
end #function (l::AddLayer)(li::Array{Layer,1})

function (l::AddLayer)(ls::Array{L,1}) where {L<:Layer}
    for li in ls
        if !in(li,l.prevLayer)
            push!(l.prevLayer, li)
        end
        if !in(l, li.nextLayers)
            push!(li.nextLayers, l)
        end
    end #for

    l.inputS = l.outputS = ls[1].outputS
    l.channels = ls[1].channels

    return l
end #function (l::AddLayer)(li::Array{Layer,1})

"""
    define input as X
"""
function (l::FCLayer)(x::Array)
    l.prevLayer = nothing
    return l
end #function (l::FCLayer)(x::Array)

function (l::Conv1D)(x::Array)
    l.prevLayer = nothing
    return l
end

function (l::Conv2D)(x::Array)
    l.prevLayer = nothing
    return l
end

function (l::Conv3D)(x::Array)
    l.prevLayer = nothing
    return l
end

function (l::Input)(X::AbstractArray{T,N}) where {T,N}
    l.A = X
    channels = size(X)[end-1]
    l.channels = channels
    return l
end #function (l::Input)(X::AbstractArray{T,N})

export l


text:
function bdh(session::Session, security::AbstractString, field::AbstractString, date_start::Date, date_end::Date;
            periodicity=nothing, # periodicitySelection option
            options=nothing, # expects key->value pairs or Dict
            verbose::Bool=false,
            timeout_milliseconds::Integer=UInt32(0),
            error_handling::ErrorHandling=Unwrap()
        )

    bdh(session, security, [field], date_start, date_end, periodicity=periodicity, options=options, verbose=verbose, timeout_milliseconds=timeout_milliseconds, error_handling=error_handling)
end

function bdh(session::Session, securities::Vector{T}, field::AbstractString, date_start::Date, date_end::Date;
            periodicity=nothing, # periodicitySelection option
            options=nothing, # expects key->value pairs or Dict
            verbose::Bool=false,
            timeout_milliseconds::Integer=UInt32(0),
            error_handling::ErrorHandling=Unwrap()
        ) where {T<:AbstractString}

    bdh(session, securities, [field], date_start, date_end, periodicity=periodicity, options=options, verbose=verbose, timeout_milliseconds=timeout_milliseconds, error_handling=error_handling)
end

"""
    bdh(session::Session, security::AbstractString, fields, date_start::Date, date_end::Date;
            periodicity=nothing, # periodicitySelection option
            options=nothing, # expects key->value pairs or Dict
            verbose::Bool=false,
            timeout_milliseconds::Integer=UInt32(0),
            error_handling::ErrorHandling=Unwrap()

Runs a query for historical data. Returns a `Vector` of named tuples.

Internally, it issues a `HistoricalDataRequest` in `//blp/refdata` service.

See also [`bds`](@ref).

# Arguments

* `fields` argument is either a single string or an array of string values.

* `options` argument expects a key->value pairs or a `Dict`.

* `periodicity` expects the string value for the `periodicitySelection` option.

# Simple query example

```julia
using BLPData, DataFrames, Dates

# opens a session
session = BLPData.Session()

# query historical data
result = BLPData.bdh(session, "IBM US Equity", ["PX_LAST", "VWAP_VOLUME"], Date(2020, 1, 2), Date(2020, 1, 30))

# format result as a `DataFrame`
df = DataFrame(result)
```

# Query with optional parameters

```julia
ticker = "PETR4 BS Equity"
field = "PX_LAST"
options = Dict(
    "periodicityAdjustment" => "CALENDAR",
    "periodicitySelection" => "DAILY",
    "currency" => "BRL",
    "pricingOption" => "PRICING_OPTION_PRICE",
    "nonTradingDayFillOption" => "ACTIVE_DAYS_ONLY",
    "nonTradingDayFillMethod" => "NIL_VALUE",
    "adjustmentFollowDPDF" => false,
    "adjustmentNormal" => true,
    "adjustmentAbnormal" => true,
    "adjustmentSplit" => true
)

# query for adjusted stock price
df = DataFrame(BLPData.bdh(session, ticker, field, Date(2019, 1, 1), Date(2019, 2, 10), options=options))
```
"""
function bdh(session::Session, security::AbstractString, fields::Vector{T}, date_start::Date, date_end::Date;
            periodicity=nothing, # periodicitySelection option
            options=nothing, # expects key->value pairs or Dict
            verbose::Bool=false,
            timeout_milliseconds::Integer=UInt32(0),
            error_handling::ErrorHandling=Unwrap()
        ) where {T<:AbstractString}

    @assert !isempty(fields) "Fields vector should not be empty."

    queue, corr_id = send_request(session, "//blp/refdata", "HistoricalDataRequest") do req
        push!(req["securities"], security)
        append!(req["fields"], fields)
        req["startDate"] = date_start
        req["endDate"] = date_end

        # optional args
        if periodicity != nothing
            req["periodicitySelection"] = periodicity
        end

        if options != nothing
            for (k, v) in options
                req[k] = v
            end
        end
    end

    fields_with_date = init_fields_with_date(fields)
    local result::Dict{String, Any} = Dict()

    for_each_response_message_element(queue, corr_id, timeout_milliseconds=timeout_milliseconds, verbose=verbose) do element
        parse_historical_data_response_into!(element, result, security, fields_with_date, error_handling)
    end

    return result[security]
end

function init_fields_with_date_into!(fields::Vector, result::Vector)
    @assert length(result) == length(fields) + 1
    result[1] = "date"
    for i in 1:length(fields)
        result[i+1] = fields[i]
    end
    return result
end

function init_fields_with_date(fields::Vector{String}) :: Vector{String}
    result = Vector{String}(undef, length(fields)+1)
    init_fields_with_date_into!(fields, result)
end

function init_fields_with_date(fields::Vector{T}) :: Vector{AbstractString} where {T<:AbstractString}
    result = Vector{AbstractString}(undef, length(fields)+1)
    init_fields_with_date_into!(fields, result)
end

"""
    bdh(session::Session, securities::Vector{T1}, fields::Vector{T2}, date_start::Date, date_end::Date;
            periodicity=nothing, # periodicitySelection option
            options=nothing, # expects key->value pairs or Dict
            verbose::Bool=false,
            timeout_milliseconds::Integer=UInt32(0),
            error_handling::ErrorHandling=Unwrap()
        ) where {T1<:AbstractString, T2<:AbstractString}

Runs a query for historical data.
Returns a `Dict` where the key is the security name and value is a `Vector` of named tuples.

Internally, BLPData will process a `ReferenceDataRequest` request for each security in parallel.
"""
function bdh(session::Session, securities::Vector{T1}, fields::Vector{T2}, date_start::Date, date_end::Date;
            periodicity=nothing, # periodicitySelection option
            options=nothing, # expects key->value pairs or Dict
            verbose::Bool=false,
            timeout_milliseconds::Integer=UInt32(0),
            error_handling::ErrorHandling=Unwrap()
        ) where {T1<:AbstractString, T2<:AbstractString}

    result = Dict()

    @sync for security in securities
        @async result[security] = bdh($session, $security, $fields, $date_start, $date_end, periodicity=$periodicity, options=$options, verbose=$verbose, timeout_milliseconds=$timeout_milliseconds, error_handling=$error_handling)
    end

    return result
end

max_stars_repo_licenses: OpenSSL MIT

max_issues_repo_licenses: OpenSSL MIT

max_forks_repo_licenses: OpenSSL MIT


text:
using StochasticProcesses
using Test

@testset "StochasticProcesses.jl" begin
    # Write your tests here.
end

max_stars_repo_name: LoganDGraham/StochasticProcesses

max_issues_repo_name: LoganDGraham/StochasticProcesses

max_forks_repo_name: LoganDGraham/StochasticProcesses


text:
function getmongokeys(fgl::FactorGraph, x::Symbol, IDs)
  cvid = -1
  # TODO -- could likely just use existing mapping already in fgl
  for id in IDs
    if Symbol(fgl.g.vertices[id[1]].label) == x
      cvid = id[2]
      break
    end
  end
  if cvid == -1
    warn("getmongokeys is not finding $(x)")
    return Dict{AbstractString, Any}(), cvid
  end
  # @show cvid
  cv = CloudGraphs.get_vertex(fgl.cg, cvid)
  if haskey(cv.properties, "mongo_keys")
    jsonstr = cv.properties["mongo_keys"]
    return JSON.parse(jsonstr), cvid
  else
    return Dict{AbstractString, Any}(), cvid
  end
end

function fetchmongorgbimg(cg::CloudGraph, key::AbstractString)
  myKeyToFind = BSONOID(key)
  findsomthing = find(cg.mongo.cgBindataCollection, ("_id" => eq(myKeyToFind)))
  myFavouriteKey = first( findsomthing );
  data = myFavouriteKey["val"]
  img = ImageMagick.readblob(data);

  # r, c = size(img)
  # imgA = zeros(r,c,3);
  # for i in 1:r, j in 1:c
  #   imgA[i,j,1] = img[i,j].r
  #   imgA[i,j,2] = img[i,j].g
  #   imgA[i,j,3] = img[i,j].b
  # end
  return img
end

function bin2arr(data::Vector{UInt8}; dtype::DataType=Float64)
  len = length(data)
  dl = sizeof(dtype)
  alen = round(Int,len/dl)

  ptrT = pointer(data)  # Ptr{UInt8}
  ptrTf = convert(Ptr{dtype}, ptrT)  # Ptr{Float32}
  arr = Vector{dtype}(alen);

  unsafe_copy!(pointer(arr),ptrTf,alen)

  return arr
end

function fetchmongodepthimg(cg::CloudGraph, key::AbstractString; dtype::DataType=Float64)
  myKeyToFind = BSONOID(key) # some valid numbers

  findsomthing = find(cg.mongo.cgBindataCollection, ("_id" => eq(myKeyToFind)))
  myFavouriteKey = first( findsomthing );
  mfkv = myFavouriteKey["val"];

  return bin2arr(mfkv, dtype=dtype)
  # len = length(mfkv)
  # dl = sizeof(dtype)
  # alen = round(Int,len/dl)
  #
  # ptrT = pointer(mfkv)  # Ptr{UInt8}
  # ptrTf = convert(Ptr{dtype}, ptrT)  # Ptr{Float32}
  # arr = Vector{dtype}(alen);
  #
  # unsafe_copy!(pointer(arr),ptrTf,alen)
  #
  # return arr
end

function cloudimshow(cg::CloudGraph,
        session::AbstractString,
        vsym::Symbol;
        descr::AbstractString = "keyframe_rgb"  )
  #
  cv = getCloudVert(cg, session, sym=vsym, bigdata=true)
  imdata = Caesar.getBigDataElement(cv, descr).data
  imshowhackpng(imdata)
end

function cloudimshow(cg::CloudGraph;
        neoid::Int=-1,
        descr::AbstractString = "keyframe_rgb"  )
  #
  cv = CloudGraphs.get_vertex(cg, neoid, true)
  imdata = Caesar.getBigDataElement(cv, descr).data
  imshowhackpng(imdata)
end

function findAllBinaryFactors(cgl::CloudGraph, session::AbstractString)
  xx = ls(cgl, session)

  slowly = Dict{Symbol, Tuple{Symbol, Symbol, Symbol}}()
  @showprogress 1 "Finding all binary edges..." for (x,va) in xx
    facts = ls(cgl, session, sym=x)
    for (fc, va2) in facts
      nodesdict = ls(cgl, session, sym=fc)
      if length(nodesdict) == 2
        # add to dictionary for later drawing
        nodes = collect(keys(nodesdict))
        if !haskey(slowly, fc) && !haselement(nodesdict[nodes[1]][3],:FACTOR) && !haselement(nodesdict[nodes[2]][3],:FACTOR)
          # fv = getVert(fgl, fgl.fIDs[fc])
          # vty = typeof(getfnctype(fv)).name.name
          slowly[fc] = (nodes[1], nodes[2], Symbol("NEEDCACHING"))
        end
      end
    end
  end

  return slowly
end

function cachepointclouds!(cache::Dict, cv::CloudVertex, ke::AbstractString, param::Dict)
  if !haskey(cache, ke)
    data = getBigDataElement(cv,ke)
    if typeof(data) == Nothing
      # warn("unable to load $(ke) from Mongo, gives data type Nothing")
      return nothing
    end
    data = data.data
    if ke == "depthframe_image"
      cache[ke] = getPointCloudFromKinect(data, param["dcamjl"], param["imshape"])
      # ri,ci = param["imshape"][1], param["imshape"][2] # TODO -- hack should be removed since depth is array and should have rows and columns stored in Mongo
      # # arrdata = data.data
      # arr = bin2arr(data, dtype=Float32) # should also store dtype for arr in Mongo
      # img = reshape(arr, ci, ri)'
      # X = reconstruct(param["dcamjl"], Array{Float64,2}(img))
      # cache[ke] = X
    elseif ke == "BSONpointcloud"
      # deserialize BSON-encoded pointcloud
      cache[ke] = getPointCloudFromBSON(data)
      # buf = IOBuffer(data)
      # st = takebuf_string(buf)
      # bb = BSONObject(st)
      # ptarr = map(x -> convert(Array, x), bb["pointcloud"])
      # cache[ke] = ptarr
    end
  end
  nothing
end

function retrievePointcloudColorInfo!(cv::CloudVertex, va::AbstractString)
  rgb = Array{Colorant,2}()
  if !hasBigDataElement(cv, va)
    warn("could not find color map in mongo, $(va)")
    return rgb
  end
  data = getBigDataElement(cv, va).data
  if va == "keyframe_rgb" || va == "keyframe_segnet"
    rgb = ImageMagick.readblob(data);
  elseif va == "BSONcolors"
    buffer = IOBuffer(data)
    str = takebuf_string(buffer)
    bb = BSONObject(str)
    # TODO -- maybe better to do: map(f, x), where f(x),
    # see http://docs.julialang.org/en/stable/manual/style-guide/#do-no-write-x-f-x
    carr = map(x -> convert(Array{UInt8}, x), bb["colors"])
    # typeof(rgb) = Array{Array{Colorant,1},1}
    rgb = carr
  end

  return rgb
end

function robotsetup(cg::CloudGraph, session::AbstractString)
  resp = fetchrobotdatafirstpose(cg, session)

  if haskey(resp, "CAMK")
    # CAMK = [[ 570.34222412; 0.0; 319.5]';
    #  [   0.0; 570.34222412; 239.5]';
    #  [   0.0; 0.0; 1.0]'];
    dcamjl = DepthCamera(resp["CAMK"])
    buildmesh!(dcamjl)
    resp["dcamjl"] = dcamjl
  end

  if haskey(resp, "bTc")
    bTc = Translation(0.0,0,0) ∘ LinearMap( CoordinateTransformations.Quat(1.0,0,0,0) )
    if resp["bTc_format"] == "xyzqwqxqyqz"
      bTc = Translation(resp["bTc"][1:3]...) ∘ LinearMap( CoordinateTransformations.Quat(resp["bTc"][4:7]...) )
    else
      warn("Unknown bTc_format, assuming identity for bTc")
    end
    resp["bTc"] = bTc
  end
  resp
end

#


text:
# export PSTH, PSTHPlot

struct PSTHResult{T}
    ROIRange::StepRangeLen
    PSTHCount::Vector{T}
    Raster::Vector{Vector{T}}
end

"""
PSTH(train::Vector{T}, marker::Vector{T}, ROI::Tuple{T,T}; binsize::T=0.1) where {T}

PSTH - peri-stimulus time histogram or post-stimulus time histogram
each bin is calcuated as \$\$\\frac{\\sum_i^N x_i}{N \\cdot \\text{binsize}}\$\$

arguments:
- train: the spike train as 1d numpy.array
- marker: the marker of stimuli
- ROI: the region of interest for each stimuli,
       define the range the PSTH.

keyword arguments:
- binsize: the size of each bin [default: 0.1]

return:
- PSTHResult{T}(ROIRange, PSTH, _seg)
    - ROIRange: the x axis
    - PSTH: the `spike/binsize` value of each PSTH bin
    - _seg: spike train for each segment, i.e. the data
            of the raster plot.

"""
function PSTH(train::Vector{T}, marker::Vector{T}, ROI::Tuple{T,T}; binsize::T=0.1) where {T}

    _seg = Vector{Vector{T}}(undef, length(marker))
    for (idx, item) in enumerate(marker)
        _seg[idx] = train[(item + ROI[1] .<= train) .& (train .< item + ROI[2])] .- item
    end

    _roi_len = Int(ceil((ROI[2] - ROI[1])/binsize))
    _roi_range = range(ROI[1], stop=ROI[1]+binsize*_roi_len, length=_roi_len+1)

    _psth = zeros(T, _roi_len)
    for item in _seg
        for i in item
            _psth[Int(floor((i - ROI[1]) / binsize)) + 1] += 1
        end
    end
    _psth = _psth ./ length(marker) ./ binsize

#     _roi_range, _psth, _seg
    PSTHResult{T}(_roi_range, _psth, _seg)
end
PSTH(sch::SpikeChannel, ech::EventChannel, ROI::Tuple{Float64, Float64}; binsize::Float64=0.1) = PSTH(sch.Value, ech.Value, ROI, binsize=binsize)

function PSTHPlot(r::PSTHResult)

    plot(legend=nothing)
    for (idx,item) in enumerate(r.Raster)
        for each in item
            plot!([each, each], [idx-0.5, idx+0.5], color=:black)
        end
    end
    _raster = plot!(ylabel="trial #")

    _step = r.ROIRange.step.hi + r.ROIRange.step.lo
    _psth = bar(r.ROIRange[1:end-1] .+ _step / 2, r.PSTHCount, legend=nothing, xlabel="ROI (sec)", ylabel="PSTH (spike/sec)")
    plot(_raster, _psth, layout=@layout([_raster; _psth]), size=(600, 800))
end


text:
using Rewrite

export unify, Unifier

struct Unifier <: AbstractDict{Term,Term}
    dict::Dict{Term,Term}
end
Unifier(ps::Pair...) = Unifier(Dict{Term,Term}(ps))
Base.length(σ::Unifier) = length(σ.dict)
Base.iterate(σ::Unifier) = iterate(σ.dict)
Base.iterate(σ::Unifier, state) = iterate(σ.dict, state)
Base.getindex(σ::Unifier, t::Term) = getindex(σ.dict, t)
function Base.setindex!(σ::Unifier, value::Term, key::Term)
    value == key && return σ

    occursin(key, value) && return nothing
    σ.dict[key] = value

    for (k, v) ∈ pairs(σ.dict)
        σ.dict[k] = σ(v)
    end

    σ
end
Base.get(σ::Unifier, t::Term, default) = get(σ.dict, t, default)

(σ::Unifier)(t::Term) = replace(t, σ)
(σ::Unifier)(xs) = map(σ, xs)

function Base.merge!(σ::Unifier, σs::Unifier...)
    for σ′ ∈ σs
        for (k, v) in σ′
            if haskey(σ, k)
                σ[k] == v || return nothing
            else
                setindex!(σ, v, k) === nothing && return nothing
            end
        end
    end
    σ
end
Base.merge(σ::Unifier, σs::Unifier...) = merge!(Unifier(), σ, σs...)

"""
    unify(t, u) -> Union{Unifier, Nothing}

Unifier terms `t` and `u`, producing a `Unifier` if the process succeeds.

# Examples
```jldoctest
julia> unify(@term(f(y)), @term(x))
Unifier with 1 entry:
  @term(x) => @term(f(y))

julia> unify(@term(x), @term(x))
Unifier with 0 entries

julia> unify(@term(x), @term(f(x)))

julia> unify(@term(f(x)), @term(g(x)))

julia> unify(@term(f(x, y)), @term(f(y, z)))
Unifier with 2 entries:
  @term(x) => @term(z)
  @term(y) => @term(z)
```
"""
unify(t::Term, u::Term) = _unify((t, u))

_unify(σ::Unifier, (x, y)::Tuple{Variable,Variable}, ms...) =
    x == y ? _unify(σ, ms...) : eliminate!(σ, (x, y), ms)
_unify(σ::Unifier, (x, t)::Tuple{Variable,Term}, ms...) = eliminate!(σ, (x, t), ms)
_unify(σ::Unifier, (t, x)::Tuple{Term,Variable}, ms...) = eliminate!(σ, (x, t), ms)
_unify(σ::Unifier, (f, g)::Tuple{Fn,Fn}, ms...) =
    f.name == g.name && length(f) == length(g) ? _unify(σ, zip(f, g)..., ms...) : nothing
_unify(σ::Unifier, (a, b)::Tuple{T,T}, ms...) where {T<:Constant} =
    get(a) == get(b) ? _unify(σ, ms...) : nothing
_unify(σ::Unifier, ms...) = nothing
_unify(σ::Unifier) = σ
_unify(ms...) = _unify(Unifier(), ms...)

function eliminate!(σ::Unifier, (x, t)::Tuple{Variable,Term}, ms)
    setindex!(σ, t, x) === nothing && return nothing
    _unify(σ, Unifier(x => t).(ms)...)
end