tBuLi / symfit

Symbolic Fitting; fitting as it should be.
http://symfit.readthedocs.org
MIT License
233 stars 17 forks source link

Questions about Fourier Series #311

Closed JohnsonW-z closed 3 years ago

JohnsonW-z commented 4 years ago

When I perform an n=100 Fourier progression calculation with an array containing 1024 elements, I get an error. How should I improve this? Doubts from a rookie. The code as follows:

`from symfit import parameters, variables, sin, cos, Fit import numpy as np import matplotlib.pyplot as plt def fourier_series(x, f, n=0): a0, cos_a = parameters(','.join(['a{}'.format(i) for i in range(0, n + 1)])) sin_b = parameters(','.join(['b{}'.format(i) for i in range(1, n + 1)])) series = a0 + sum(ai cos(i f x) + bi sin(i f * x) for i, (ai, bi) in enumerate(zip(cos_a, sin_b), start=1)) return series x, y = variables('x, y') w, = parameters('w') model_dict = {y: fourier_series(x, f=w, n=3)} print(model_dict)

xdata = np.array([0.0076029,0.007621867,0.007642907,0.007665746,0.007690046,0.007715406,0.007741381,0.007767494,0.007793251,0.007818159,0.007841746,0.007863574,0.007883263,0.007900508,0.007915093,0.007926906,0.007935945,0.007942316,0.007946232,0.007947997,0.007947993,0.00794666,0.00794448,0.007941954,0.007939586,0.007937866,0.007937249,0.007938146,0.007940904,0.007945798,0.007953026,0.0079627,0.007974856,0.007989459,0.008006412,0.008025571,0.008046759,0.008069774,0.008094408,0.008120447,0.008147684,0.008175921,0.008204973,0.008234669,0.008264855,0.008295392,0.008326155,0.008357034,0.008387931,0.008418755,0.00844943,0.008479881,0.008510045,0.008539864,0.008569288,0.008598276,0.008626797,0.008654834,0.008682384,0.008709459,0.008736089,0.008762319,0.008788211,0.00881384,0.008839292,0.008864661,0.008890047,0.00891555,0.008941264,0.008967278,0.008993667,0.009020486,0.009047769,0.009075523,0.009103726,0.009132328,0.009161248,0.009190381,0.009219599,0.00924876,0.00927771,0.009306291,0.009334347,0.009361734,0.009388322,0.009414005,0.009438709,0.009462395,0.009485065,0.009506767,0.009527588,0.00954766,0.00956715,0.009586254,0.009605189,0.009624182,0.009643461,0.009663239,0.0096837,0.009704985,0.009727173,0.009750271,0.009774192,0.009798756,0.009823678,0.009848573,0.009872963,0.009896288,0.009917926,0.009937211,0.009953454,0.009965976,0.009974128,0.009977332,0.009975108,0.009967109,0.009953154,0.009933247,0.009907592,0.009876597,0.009840861,0.009801148,0.00975836,0.009713488,0.009667569,0.009621625,0.00957661,0.009533334,0.009492398,0.009454116,0.009418444,0.00938492,0.009352613,0.009320106,0.009285489,0.009246395,0.009200051,0.009143352,0.009072962,0.00898542,0.008877282,0.008745266,0.008586431,0.008398357,0.008179336,0.007928546,0.007646195,0.00733361,0.006993266,0.006628743,0.006244617,0.005846298,0.005439828,0.005031649,0.004628346,0.004236376,0.003861774,0.003509841,0.003184835,0.002889669,0.002625651,0.002392294,0.002187224,0.002006185,0.001843169,0.001690628,0.001539764,0.001380862,0.001203657,0.000997722,0.000752868,0.000459551,0.000109281,-0.000304994,-0.000788562,-0.001344722,-0.001974707,-0.002677782,-0.003451535,-0.004292348,-0.005196016,-0.006158455,-0.007176458,-0.00824842,-0.009374973,-0.010559462,-0.011808195,-0.013130388,-0.014537807,-0.016044087,-0.017663844,-0.01941172,-0.021301539,-0.023345794,-0.025555569,-0.027940953,-0.030511825,-0.033278814,-0.036254139,-0.039452078,-0.042888923,-0.046582469,-0.050551323,-0.054814503,-0.059391897,-0.064306067,-0.069585631,-0.075270033,-0.081415014,-0.088097695,-0.095419936,-0.103508714,-0.112512674,-0.122594652,-0.13392073,-0.146647029,-0.160905817,-0.17679258,-0.194355386,-0.213587472,-0.234423514,-0.256739698,-0.280357508,-0.305050975,-0.330556951,-0.356587737,-0.382845052,-0.409034194,-0.434877148,-0.460123611,-0.484559234,-0.508010798,-0.530348366,-0.551484725,-0.571372491,-0.589999333,-0.607381811,-0.623558391,-0.638582266,-0.652514634,-0.665418982,-0.677356771,-0.688384667,-0.698553229,-0.707906753,-0.716483887,-0.724318584,-0.731441063,-0.737878537,-0.743655631,-0.748794494,-0.753314711,-0.757233118,-0.760563637,-0.763317202,-0.765501812,-0.76712273,-0.768182805,-0.768682879,-0.768622249,-0.767999123,-0.766811016,-0.765055034,-0.762727987,-0.759826293,-0.756345671,-0.75228067,-0.747624085,-0.742366375,-0.736495158,-0.729994857,-0.722846509,-0.715027699,-0.706512563,-0.697271748,-0.687272271,-0.676477261,-0.664845651,-0.652332013,-0.638886783,-0.624457195,-0.608989191,-0.592430523,-0.574735064,-0.55586821,-0.535813047,-0.514576828,-0.492197226,-0.468747805,-0.444342139,-0.419136021,-0.393327191,-0.367152061,-0.340879097,-0.314798861,-0.289211168,-0.264410387,-0.240670245,-0.218229691,-0.197281182,-0.177962405,-0.160351913,-0.144468771,-0.130275964,-0.117687174,-0.106576384,-0.09678957,-0.088157507,-0.080508505,-0.07367977,-0.06752622,-0.061926005,-0.056782472,-0.052022968,-0.047595333,-0.043463188,-0.039601113,-0.035990562,-0.03261697,-0.029468095,-0.026533326,-0.023803505,-0.021270784,-0.018928206,-0.016768875,-0.014784843,-0.012966001,-0.01129935,-0.00976896,-0.008356769,-0.007044153,-0.005813945,-0.00465244,-0.003550886,-0.002506041,-0.001519622,-0.00059673,0.000256402,0.001034819,0.001736859,0.002365344,0.002927715,0.003435155,0.003900964,0.004338618,0.004759924,0.005173672,0.005584943,0.005995138,0.006402563,0.006803379,0.007192649,0.007565301,0.007916855,0.008243881,0.008544191,0.008816822,0.009061892,0.009280381,0.009473893,0.009644438,0.009794242,0.009925589,0.010040709,0.010141695,0.010230446,0.01030864,0.010377728,0.010438941,0.010493307,0.010541681,0.01058477,0.010623162,0.010657351,0.010687753,0.010714723,0.010738563,0.010759528,0.01077783,0.010793638,0.010807078,0.010818239,0.01082717,0.01083389,0.010838387,0.010840629,0.010840565,0.010838138,0.010833285,0.010825946,0.010816071,0.010803625,0.010788592,0.010770985,0.010750843,0.01072824,0.010703278,0.01067609,0.010646835,0.010615691,0.010582855,0.010548531,0.010512928,0.010476256,0.010438714,0.010400492,0.010361761,0.010322671,0.010283347,0.010243888,0.010204364,0.01016482,0.010125275,0.010085731,0.010046174,0.01000658,0.009966918,0.009927158,0.00988727,0.00984723,0.009807021,0.009766638,0.009726085,0.00968538,0.009644555,0.009603655,0.009562739,0.009521878,0.009481156,0.009440663,0.009400501,0.009360774,0.009321587,0.009283043,0.009245238,0.009208256,0.009172164,0.009137009,0.009102812,0.009069574,0.009037266,0.00900584,0.00897523,0.008945353,0.008916117,0.008887424,0.008859174,0.008831268,0.008803611,0.008776113,0.00874869,0.008721267,0.008693776,0.008666158,0.008638362,0.008610343,0.008582063,0.008553494,0.00852461,0.008495395,0.008465841,0.00843595,0.00840574,0.008375246,0.008344527,0.008313673,0.008282804,0.008252079,0.0082217,0.008191908,0.008162985,0.008135247,0.008109043,0.008084744,0.008062729,0.008043376,0.008027041,0.008014038,0.008004617,0.007998936,0.007997043,0.00799885,0.008004123,0.008012475,0.008023371,0.008036137,0.008049983,0.00806403,0.008077336,0.008088934,0.008097865,0.008103222,0.008104188,0.008100082,0.008090401,0.008074856,0.008053403,0.008026259,0.007993898,0.007957046,0.007916644,0.007873818,0.007829829,0.00778603,0.007743813,0.007704558,0.007669581,0.007640082,0.007617091,0.007601426,0.007593651,0.007594053,0.007602632,0.007619113,0.007642964,0.007673443,0.007709636,0.007750508,0.007794949,0.007841824,0.007890007,0.007938423,0.007986078,0.008032091,0.008075712,0.008116341,0.008153535,0.008187003,0.008216601,0.008242305,0.008264198,0.008282438,0.008297238,0.008308842,0.00831751,0.008323504,0.008327075,0.008328463,0.008327885,0.008325541,0.008321607,0.008316242,0.008309582,0.008301749,0.008292847,0.008282967,0.008272186,0.008260575,0.008248192,0.008235094,0.008221331,0.008206949,0.008191996,0.008176516,0.008160555,0.008144158,0.008127369,0.008110235,0.008092799,0.008075106,0.008057199,0.008039121,0.008020914,0.00800262,0.007984281,0.007965938,0.007947636,0.007929415,0.007911322,0.0078934,0.007875694,0.00785825,0.007841112,0.007824321,0.007807918,0.007791937,0.007776403,0.007761339,0.007746755,0.007732653,0.007719028,0.007705866,0.007693147,0.007680846,0.007668937,0.00765739,0.007646178,0.007635273,0.007624648,0.007614281,0.007604151,0.00759424,0.007584532,0.007575016,0.007565682,0.007556521,0.007547527,0.007538699,0.007530032,0.007521525,0.007513179,0.007504992,0.007496963,0.007489091,0.007481374,0.007473807,0.007466387,0.007459107,0.00745196,0.00744494,0.007438038,0.007431248,0.007424563,0.007417976,0.007411481,0.007405074,0.00739875,0.007392508,0.007386344,0.007380258,0.00737425,0.00736832,0.00736247,0.007356703,0.00735102,0.007345425,0.007339921,0.007334509,0.007329193,0.007323972,0.007318846,0.007313816,0.007308877,0.007304026,0.00729926,0.007294572,0.007289956,0.007285407,0.007280918,0.007276483,0.007272097,0.007267754,0.007263449,0.007259177,0.007254934,0.007250716,0.007246519,0.007242341,0.007238179,0.007234029,0.007229889,0.007225757,0.00722163,0.007217508,0.007213388,0.007209269,0.007205149,0.007201028,0.007196903,0.007192775,0.007188643,0.007184505,0.007180361,0.00717621,0.007172052,0.007167885,0.007163709,0.007159523,0.007155325,0.007151114,0.007146889,0.007142649,0.00713839,0.007134112,0.007129812,0.007125489,0.007121139,0.007116761,0.007112353,0.007107914,0.007103442,0.007098937,0.0070944,0.007089832,0.007085235,0.007080613,0.00707597,0.007071311,0.007066642,0.007061969,0.007057301,0.007052643,0.007048005,0.007043395,0.007038822,0.007034296,0.007029827,0.007025427,0.007021107,0.00701688,0.00701276,0.007008759,0.007004893,0.007001175,0.006997618,0.006994235,0.006991036,0.006988029,0.00698522,0.006982609,0.006980195,0.006977971,0.006975929,0.006974056,0.006972337,0.006970757,0.006969299,0.006967945,0.006966677,0.006965478,0.006964328,0.006963212,0.006962111,0.006961008,0.006959887,0.006958728,0.006957516,0.006956234,0.006954863,0.006953388,0.006951794,0.006950065,0.00694819,0.00694616,0.006943967,0.00694161,0.00693909,0.006936412,0.006933585,0.006930622,0.006927537,0.006924347,0.006921071,0.00691773,0.006914343,0.006910933,0.00690752,0.006904126,0.006900775,0.00689749,0.006894293,0.006891212,0.006888269,0.006885492,0.006882907,0.006880539,0.006878415,0.006876556,0.006874986,0.006873719,0.006872769,0.006872141,0.006871837,0.006871849,0.006872166,0.006872771,0.006873641,0.006874751,0.006876075,0.006877583,0.006879247,0.006881036,0.006882921,0.006884873,0.006886864,0.006888865,0.006890847,0.006892782,0.006894641,0.006896396,0.006898019,0.006899482,0.006900759,0.006901824,0.006902656,0.006903237,0.006903552,0.006903593,0.00690336,0.006902855,0.006902089,0.006901077,0.00689984,0.006898402,0.006896789,0.006895029,0.006893153,0.00689119,0.00688917,0.006887124,0.006885084,0.00688308,0.006881144,0.006879308,0.006877605,0.006876068,0.006874731,0.006873626,0.006872787,0.006872245,0.006872027,0.006872157,0.006872654,0.006873531,0.006874792,0.006876435,0.006878449,0.006880818,0.00688352,0.006886528,0.006889812,0.006893339,0.006897076,0.006900991,0.00690505,0.00690922,0.006913468,0.006917765,0.006922077,0.006926376,0.006930629,0.006934807,0.006938879,0.006942816,0.00694659,0.006950174,0.006953543,0.006956675,0.006959554,0.006962167,0.006964509,0.006966579,0.006968383,0.006969932,0.006971244,0.006972338,0.006973238,0.006973969,0.006974557,0.00697503,0.006975414,0.006975738,0.006976028,0.006976311,0.006976614,0.006976965,0.00697739,0.006977918,0.006978575,0.006979389,0.006980388,0.006981596,0.006983039,0.006984736,0.006986708,0.006988965,0.006991516,0.006994365,0.006997506,0.007000931,0.007004628,0.007008579,0.007012764,0.007017162,0.007021751,0.007026508,0.007031411,0.007036438,0.007041569,0.007046783,0.007052062,0.007057386,0.007062739,0.007068103,0.007073462,0.0070788,0.007084101,0.007089352,0.007094539,0.007099649,0.007104674,0.007109605,0.007114437,0.007119166,0.007123793,0.007128319,0.00713275,0.007137091,0.00714135,0.007145535,0.007149656,0.007153722,0.007157743,0.007161728,0.007165685,0.007169624,0.007173554,0.007177482,0.007181417,0.007185367,0.007189339,0.007193343,0.007197384,0.007201472,0.007205613,0.007209816,0.007214086,0.007218431,0.007222856,0.007227368,0.007231971,0.007236671,0.007241472,0.00724638,0.007251399,0.007256535,0.007261796,0.007267188,0.007272722,0.007278408,0.007284259,0.00729029,0.007296519,0.007302964,0.007309648,0.007316594,0.00732383,0.007331383,0.007339283,0.007347559,0.007356239,0.007365351,0.007374916,0.007384951,0.007395467,0.007406466,0.00741794,0.007429875,0.007442248,0.007455029,0.007468182,0.007481669,0.007495447,0.007509473,0.007523704,0.007538094,0.007552602,0.007567182,0.007581792,0.00759639,0.007610932,0.007625376,0.007639678,0.007653795,0.007667684,0.007681303,0.007694612,0.007707573,0.007720152,0.007732322,0.007744062,0.007755358,0.007766207,0.007776612,0.007786585,0.007796144,0.007805312,0.007814118,0.00782259,0.007830758,0.007838653,0.007846305,0.007853741,0.007860987,0.007868068,0.007875005,0.007881819,0.007888529,0.007895152,0.007901703,0.007908198,0.00791465,0.007921071,0.007927472,0.007933865,0.007940259,0.007946662,0.007953084,0.007959531,0.00796601,0.007972529,0.007979092,0.007985705,0.007992376,0.007999109,0.008005913,0.008012794,0.00801976,0.008026822,0.008033989,0.008041274,0.00804869,0.008056252,0.008063976,0.008071881,0.008079985,0.008088308,0.008096873,0.0081057,0.008114809,0.00812422,0.008133948,0.008144005,0.0081544 ]) ydata = np.linspace(0, 2*np.pi,1024)

fit = Fit(model_dict, x=xdata, y=ydata) fit_result = fit.execute() print(fit_result) plt.plot(xdata, ydata) plt.plot(xdata, fit.model(x=xdata, **fit_result.params).y, ls=':') plt.xlabel('x') plt.ylabel('y') plt.show()`

pckroon commented 4 years ago

Hey,

first, could you edit your post to use triple ` instead of single so the code gets formatted nicely? Secondly, what error do you get?

tBuLi commented 4 years ago

I wouldn't be surprised if it is some kind of "about to explode" error ;). With a hundred coefficients this becomes a huge model, this example was never intended to be used that way because it's terribly inefficient for serious n. Perhpas you can comment on what it is you are trying to do? I think using a discrete FFT might be a better option?

JohnsonW-z commented 4 years ago

I wouldn't be surprised if it is some kind of "about to explode" error ;). With a hundred coefficients this becomes a huge model, this example was never intended to be used that way because it's terribly inefficient for serious n. Perhpas you can comment on what it is you are trying to do? I think using a discrete FFT might be a better option?

I expect to get the value of each 'a' and 'b' (n from 1 to 100) by using the Fourier transform.But I could never find a way. I'm too weak. sad QAQ

JohnsonW-z commented 4 years ago

Hey,

first, could you edit your post to use triple ` instead of single so the code gets formatted nicely? Secondly, what error do you get?

This is my first post and don't know how to paste the code. clicked the "insert code" button .I will pay attention next time. '''MemoryError: Unable to allocate array with shape (202, 202, 1024) and data type `float64' '''

pckroon commented 4 years ago

No worries, there's a first time for everything :) Welcome aboard.

You can even edit your original post to adapt the formatting (and check out the preview button). The single backticks are for inline code: fit = Fit(...). Triple backticks are for code blocks. See also 1 and 2

As for the actual error, in the future also post the full traceback you get, rather than just the final line. The traceback might tell us where things are going wrong (and by extension why). In this case though, it's very likely that your model is simply too large (too many parameters) for the hardware you have available. That said, you can try fitting with a derivative-free minimizer such as Nelder Mead:

from symfit.core.minimizers import NelderMead
...
fit = Fit(..., minimizer=NelderMead)

The better option is still probably to do an actual FFT, using either fft, rfft or hfft from e.g. https://numpy.org/doc/stable/reference/routines.fft.html. Also take a look at the helper functions there.

tBuLi commented 3 years ago

By the time you have this many parameters you should definitely switch to a fourier transform. You might also like my package Deconvoluted, which I made because I think the API of Fourier transforms in numpy/scipy is a bit... well, convoluted. Checkout the 1D example here, might be close to what you need: https://deconvoluted.readthedocs.io/en/latest/examples/fourier_1d.html